# Random Finite Set Based Data Assimilation for Dynamic Data Driven Simulation of Maritime Pirate Activity

Authors: Peng Wang; Ge Li; Rusheng Ju; Yong Peng
Journal: Mathematical Problems in Engineering
 (2017)

---

## Abstract

Maritime piracy is posing a genuine threat to maritime transport. The main purpose of simulation is to predict the behaviors of many actual systems, and it has been successfully applied in many fields. But the application of simulation in the maritime domain is still scarce. The rapid development of network and measurement technologies brings about higher accuracy and better availability of online measurements. This makes the simulation paradigm named as dynamic data driven simulation increasingly popular. It can assimilate the online measurements into the running simulation models and ensure much more accurate prediction of the complex systems under study. In this paper, we study how to utilize the online measurements in the agent based simulation of the maritime pirate activity. A new random finite set based data assimilation algorithm is proposed to overcome the limitations of the conventional vectors based data assimilation algorithms. The random finite set based general data model, measurement model, and simulation model are introduced to support the proposed algorithm. The details of the proposed algorithm are presented in the context of agent based simulation of maritime pirate activity. Two groups of experiments are used to practically prove the effectiveness and superiority of the proposed algorithm.


---

## Body

## 1. Introduction

In modern society, global maritime transport system is a critical part of the world’s transportation system. According to information, nearly ninety percent of internationally traded goods are transported by sea [1]. Since 2010, insurance rate has increased more than 10 times for ships transiting the areas where the maritime piracy is rampant. In 2016, the insured losses caused by maritime piracy worldwide amounted to 2.46 billion US dollars and are rising continually. The number of pirate attacks against ships worldwide in 2016 was 191. The increasingly rampant maritime pirate activities have become a great threat to the international maritime transport system throughout the world.Since entering the new century, many researchers have studied how to fight against the dramatically increased maritime piracy [2]. Many of them offer insight into how the problem may be addressed through global political, security studies, international relations, or some other approaches. These approaches proposed by them target the restoration of order and law in weak states, but they offer only the most general intelligence while state fragility still remains [3–5]. For example, many approaches only aim to tell where attacks are most likely to happen without giving advice on how to prevent or stop them.The main purpose of simulation is to predict the behaviors of many actual systems, such as natural disasters, wildfires, and terrorist attacks. The simulation systems have been used to support policy design and operational management in the ground transportation [6] and air transportation field [7] for a very long time. But the applications of simulation systems in the maritime transportation field are scarce.The agent based dynamic data driven simulation (DDDS) of maritime pirate activity focuses on the accurate prediction of marine pirates in the real world and takes full advantage of the real-time sensors. The primary purpose of the simulation system is to accurately assess the potential outcomes of actions taken by decision-makers (including ship commanders and operators) who interact with the system. It also can be used to find the real-time optimal penetration routes through designated piracy zones by assimilating the real-time measurements from all kinds of sensors. Data assimilation techniques, which are the foundation of DDDS, can fuse real-time measurements with simulated states to best determine the states of the studied system. Since actions are chosen according to some policy based on the current simulation states of the system, simulation is deservedly used to determine the utility of these actions. The facts that make the data assimilation technology necessary for DDDS of maritime pirate activity are as follows.Highly Dynamic and Changing Real Operational Context. Future maritime mission planning environments are expected to be complex and distributed due to rapid advancement in cyber-physical systems and the ubiquitous use of intelligent heterogeneous assets. The maritime environmental context is always changing dynamically. The simulation results based on current data may no longer be accurate as additional data is received. Unforeseen pirates can appear in the designated piracy zone and critical rescue equipment can be damaged, thereby making the current simulation results useless. In order to adapt to the dynamic and changing real operational context, we need to use the data assimilation technology to assimilate the changing environmental information into the simulation system.Imprecise Simulation Data and Models. The accuracy of simulation is determined by many factors, such as data used for initialization and models for describing the real system. Taking the simulation of maritime pirate activity as example, the accuracy of simulation depends on the weather data, merchant vessel data, and piracy vessel data in designated piracy zones. Due to current technical constraints, it is impossible to obtain all the required data without errors. Besides data errors, the simulation model of pirates also brings errors. Being subject to these constraints, the simulation results of the simulation models will certainly be inaccurate. In order to overcome the adverse impacts of the imprecise input data and models, we need the data assimilation technology.The differences between the simulation system and actual system will continue to expand, if we do not assimilate measurements from the actual system and dynamically adjust the simulation results. We propose a random finite set (RFS) based recursive algorithm for online assimilation of measurements in the presence of noises, data association uncertainty, detection uncertainty, and clutters. We model the simulation states and measurements as RFS and use the RFS based Bayesian inference recursion to assimilate the real-time measurements in time.The structure of the paper is organized as follows. Section2 introduces two related works and highlights the gaps of the literature covered by the current algorithms. Section 3 indicates the main RFS based models (including general data model, measurement model, and simulation model) which are the foundation of the proposed data assimilation algorithm. Section 4 introduces the agent based models and the combination of data assimilation and agent. Section 5 presents the detailed design and implementation of the RFS based data assimilation algorithm. Section 6 mentions the experiment design and the results obtained from the experiments. Section 7 presents the conclusions.


## 2. Related Works

### 2.1. Maritime Piracy Simulation

The applications of simulation systems in the maritime transportation domain are especially scarce. The majority of existing works only focus on traffic in national, coastal waters and ports. They tend to use high level equation based models for modeling maritime piracy [8]. A little work has been done to optimize ship routing in the maritime domain, but it did not take the security aspect into consideration [9].Some researchers attempt to apply mathematical modeling and computer simulation to study the maritime piracy in recent years. In [10], the authors use the discrete event simulator PANOPEA to model the maritime pirates around the Gulf of Aden. They have evaluated the effectiveness of some command and control models. The limitation is that they only consider the application for a specified scenario, and the method and models are not applicable in other scenarios.In [11], the author employs intelligence reports, meteorological forecasts, and historical pirate incidents to describe Next-generation Piracy Performance Surface (PPSN) model and mainly focuses on predicting the areas around the Horn of Africa. In [12], the authors introduce the Pirate Attack Risk Surface (PARS) model by using the probabilistic models to further improve the PPSN model. Both these two models are numerical models. Since they do not directly model behaviors or interactions of individual ships, their ability for what-if type of simulation and prediction is limited.In [13], the authors use the agent based modeling framework named as MANA (Map-Aware, Non-Uniform Automata) to find the factors that are important for the escort through the Gulf of Aden. MANA is used for defense related modeling and simulation by specifying the properties of agents. They model the escorting scenario on a tactical level by using MANA. In [14], the authors use the MANA framework to defend vessels against pirate activities by employing nonlegal deterrents. They model the single encounter behaviors in detail, and have not taken the whole system into consideration. The limitation is that it has not taken the measurements of real situation into consideration, which makes it unpractical in applications.Jakob et al. are the first to use agent based simulation to study the international maritime transport security [15]. In [16], they describe and model the transportation system from a macro perspective. They describe the pirates and piracy countermeasures by using the agent based models. And they drive the agent model by inputting the offline data stored in the database. They conduct the simulation with thousands of agents, and each agent is related to a vessel. They can capture the dynamic of the real transportation system and evaluate the effectiveness of piracy countermeasures. But they only focus on the whole maritime transportation system from a macro perspective but fail to support decision for the ship operator in a tactical level. The data used to drive the simulation are offline data that are stored in the database. They use the stored data to drive the simulation instead of real-time measurements, so the simulation cannot be used to support real-time decision-making.The agent based DDDS of maritime pirate activity studies the problem of the merchant vessel agent trying to cross or penetrate a piracy-infested area patrolled by one or more pirate vessel agents by assimilating the online measurements into the simulation models. The merchant vessel agent in the simulation can choose its self-relief measures. The simulation system can be used to help minimize the chance of hostile encounter or being hijacked and to avoid detection and interception. So the simulation system can be used by the merchant vessel operators or commanders to penetrate a piracy-infested area. This simulation also can be applied in reality, such as logistics in illegal border crossing, insecure regions, and/or smuggling interdiction. The schematic diagram of the simulation system based decision support system for maritime counter-piracy is as shown in Figure1.Figure 1
Agent and simulation based decision support system for maritime counter-piracy.Many real-world problems would benefit from the ability to assimilate measurements and predict the state of complex systems in real time. Since many complex systems, such as modeling counter-piracy operations, can take on very different states at critical decision points in the simulation, it is necessary for these simulation systems to support branching through either the generation of separate Monte Carlo replications or the use of multiple hypotheses branching techniques that support multiple hypotheses within a single execution. Multiple hypotheses branching can also be used to explore what-if possibilities or alternative plans. The traditional estimation and prediction control theory is combined with optimistic simulation techniques to support the data assimilation of complex systems. Real-time simulation based decision support system for maritime counter-piracy projects the future through simulation. It uses real-time measurements to calibrate the current state and requires multiple hypothesis what-if branching capabilities to simultaneously explore the effectiveness of various decision options.Figure2 is the schematic diagram of data assimilation for agent based DDDS of maritime pirate activity. In the simulation based decision support system, the newly assimilated simulation states of the real system at time tk are obtained by optimally combining noisy and cluttered measurements with the previously predicted state. Then, the new state estimated is extrapolated to time tk+1 to form the next prediction of the system. The prediction/update cycle repeats as time evolves and new measurements become available. This iterative process is shown in Figure 2, where X represents the abstract state of the system. Accurate assimilation must balance the uncertainty of available measurements with unknown or unpredictable behaviors associated with the model.Figure 2
The relationship between data assimilation and simulation.


### 2.2. DDDS and Data Assimilation

Until now, the simulation models have been using the static information as initial settings to simulate the real systems in the conventional simulation systems. In the dynamic maritime pirate activity modeling, since the static initial parameters fail to capture the real-time environmental changes timely, the simulation states are always quite different from the actual ones of the real system. DDDS has a golden opportunity for development, for the online measurements can reflect the latest updated states of the real system. It means dynamically assimilating online measurements when executing simulations. This idea has been applied in some research areas, such as traffic management, disaster forecasting, environmental science, and wildfire [17, 18]. The simulation system in DDDS paradigm is continually adjusted by assimilating real-time measurements for much more accurate prediction of the studied systems [19]. The rapid development of network and measurement technologies contributes to the higher accuracy and better availability of online measurements, which has facilitated the development of DDDS.DDDS highly depends on data assimilation technology which is in charge of assimilating online measurements into a running simulation system. Data assimilation has the potential to greatly improve the accuracy of simulation systems. As the simulation models become increasingly complicated, if we still make little usage of online measurements from the studied system, traditional simulations will be largely decoupled from real systems. In data assimilation, the measurements are assimilated into the running simulation models to generate more accurate simulation results [20]. The data assimilation algorithms aim to minimize the differences between the simulated states and the real world. There are two main types of data assimilation algorithms: one is the constant statistical algorithms which include 3DVAR [21] and 4DVAR [22] and the other one is the adaptive statistical algorithms, such as Kalman filter [23], ensemble Kalman filter (EnKF) [24], particle filter (PF) [25], and hierarchical Bayesian method [26]. The data assimilation algorithms that have been successfully used in simulation applications include 3DVAR, 4DVAR, Kalman filter, and PF based algorithms.The basic principle of 3DVAR and 4DVAR is to avoid the explicit calculation of the gain matrix and make its inversion by using a minimization procedure of the cost function. They are cost function based and use a suited descent algorithm to iteratively evaluate the cost function. The limitation of variational methods is that they use the minimization of a cost function and they do not provide any estimate of predictive uncertainty. Kalman filter and its variance [27, 28] are important and conventional Bayes estimation techniques used in data assimilation. Kalman filter based data assimilation algorithms are optimal for the assimilation of linear Gaussian systems. The requirements of using Kalman filter based data assimilation algorithms are that the studied systems should be modeled in linear equations and the distributions of processing noises and observation noises are Gaussian. The constraints greatly restrict the applicability of Kalman filter based algorithms in data assimilation for systems that are nonlinear and non-Gaussian. In this paper, we extend the Kalman filter by using RFS based models and unscented transformation.These conventional data assimilation algorithms fail to support data assimilation in the simulation field because of the nonlinear, non-Gaussian properties of a lot of simulation models. The sequential Monte Carlo method is also named as particle filter (PF). PF based data assimilation has been used for some simulation applications including traffic and wildfire. Particle filter is the exact method for the nonlinear and non-Gaussian applications. It has been the most widely used data assimilation algorithm in simulation applications in the recent years [19]. PF based data assimilation in traffic simulation is presented by Xie et al. in [29] and Wu et al. in [30]. PF based data assimilation in wildfire simulation is presented by Xue et al. in [31, 32]. They all use the same nonparametric statistic inference method based on particle filter, because the particle filter based data assimilation algorithms need no assumptions about the distribution or linearity of the studied simulation systems.PF based data assimilation algorithms are widely used, because they are the exact methods for the nonlinear and non-Gaussian systems. But they have the following essential disadvantages:(i)
They are based on the assumption that the studied systems are single dynamic systems that are permanently active. What is more, they use the presumed physical models of individual agents, and the transition of agents from one model to another is less common.(ii)
They are based on the assumption that there is exactly one measurement at each time and the detection is perfect with no false detections or misdetections. In other words, they assume that the sensors never fail to make a measurement of the state and never pick up false measurements.(iii)
They are based on the assumption that the measurements are precise and the measurement models are always precise as well. They always model the measurements as points in the measurement space, and the number and ordering of measurements are previously designated.The standard PF based data assimilation algorithms cannot be used for data assimilation of the dynamic systems that are switching on and off randomly. Switching is quite common in maritime pirate activity, as it did in the case where moving pirate vessels appear and vanish in the designated areas. The restriction in which the studied system must be a single dynamic system limits the applicability of the PF based algorithms. In the maritime pirate activity, multiple dynamic pirate vessels can appear and vanish during the measurement intervals. The measurement models used in PF based data assimilation algorithm are quite unreasonable, since there are always false detections and misdetections in the real world, and the measurement functions are impossible to be precisely known. The PF based algorithm cannot jointly estimate the number of pirates and the states of every pirate either.For the agent based simulation of the maritime pirate activity, we employ RFS to assimilate online measurements into the simulation model. A major task for the data assimilation is to dynamically update the simulation states of the real pirate vessels and then provide the updated states to the simulation model for follow-on simulations. The proposed RFS based data assimilation algorithm can overcome the limitations of the standard PF based data assimilation algorithms effectively.The RFS theory proposed by Goodman et al. as finite set statistics (FISST) is a recently developed theory that unifies much of information fusion under a single Bayesian paradigm [33, 34]. FISST is the first systematic treatment of multisensor multitarget filtering as part of a unified framework for data fusion using random set theory [35, 36]. FISST is a means of directly extending the familiar methodology of single-sensor, single-target Bayes statistics to multisensor multitarget problems. FISST gives some applicable mathematical tools that enable direct application of Bayesian inference. FISST provides the concept of integration/density for set-valued variable. More details can be found in [37].


## 2.1. Maritime Piracy Simulation

The applications of simulation systems in the maritime transportation domain are especially scarce. The majority of existing works only focus on traffic in national, coastal waters and ports. They tend to use high level equation based models for modeling maritime piracy [8]. A little work has been done to optimize ship routing in the maritime domain, but it did not take the security aspect into consideration [9].Some researchers attempt to apply mathematical modeling and computer simulation to study the maritime piracy in recent years. In [10], the authors use the discrete event simulator PANOPEA to model the maritime pirates around the Gulf of Aden. They have evaluated the effectiveness of some command and control models. The limitation is that they only consider the application for a specified scenario, and the method and models are not applicable in other scenarios.In [11], the author employs intelligence reports, meteorological forecasts, and historical pirate incidents to describe Next-generation Piracy Performance Surface (PPSN) model and mainly focuses on predicting the areas around the Horn of Africa. In [12], the authors introduce the Pirate Attack Risk Surface (PARS) model by using the probabilistic models to further improve the PPSN model. Both these two models are numerical models. Since they do not directly model behaviors or interactions of individual ships, their ability for what-if type of simulation and prediction is limited.In [13], the authors use the agent based modeling framework named as MANA (Map-Aware, Non-Uniform Automata) to find the factors that are important for the escort through the Gulf of Aden. MANA is used for defense related modeling and simulation by specifying the properties of agents. They model the escorting scenario on a tactical level by using MANA. In [14], the authors use the MANA framework to defend vessels against pirate activities by employing nonlegal deterrents. They model the single encounter behaviors in detail, and have not taken the whole system into consideration. The limitation is that it has not taken the measurements of real situation into consideration, which makes it unpractical in applications.Jakob et al. are the first to use agent based simulation to study the international maritime transport security [15]. In [16], they describe and model the transportation system from a macro perspective. They describe the pirates and piracy countermeasures by using the agent based models. And they drive the agent model by inputting the offline data stored in the database. They conduct the simulation with thousands of agents, and each agent is related to a vessel. They can capture the dynamic of the real transportation system and evaluate the effectiveness of piracy countermeasures. But they only focus on the whole maritime transportation system from a macro perspective but fail to support decision for the ship operator in a tactical level. The data used to drive the simulation are offline data that are stored in the database. They use the stored data to drive the simulation instead of real-time measurements, so the simulation cannot be used to support real-time decision-making.The agent based DDDS of maritime pirate activity studies the problem of the merchant vessel agent trying to cross or penetrate a piracy-infested area patrolled by one or more pirate vessel agents by assimilating the online measurements into the simulation models. The merchant vessel agent in the simulation can choose its self-relief measures. The simulation system can be used to help minimize the chance of hostile encounter or being hijacked and to avoid detection and interception. So the simulation system can be used by the merchant vessel operators or commanders to penetrate a piracy-infested area. This simulation also can be applied in reality, such as logistics in illegal border crossing, insecure regions, and/or smuggling interdiction. The schematic diagram of the simulation system based decision support system for maritime counter-piracy is as shown in Figure1.Figure 1
Agent and simulation based decision support system for maritime counter-piracy.Many real-world problems would benefit from the ability to assimilate measurements and predict the state of complex systems in real time. Since many complex systems, such as modeling counter-piracy operations, can take on very different states at critical decision points in the simulation, it is necessary for these simulation systems to support branching through either the generation of separate Monte Carlo replications or the use of multiple hypotheses branching techniques that support multiple hypotheses within a single execution. Multiple hypotheses branching can also be used to explore what-if possibilities or alternative plans. The traditional estimation and prediction control theory is combined with optimistic simulation techniques to support the data assimilation of complex systems. Real-time simulation based decision support system for maritime counter-piracy projects the future through simulation. It uses real-time measurements to calibrate the current state and requires multiple hypothesis what-if branching capabilities to simultaneously explore the effectiveness of various decision options.Figure2 is the schematic diagram of data assimilation for agent based DDDS of maritime pirate activity. In the simulation based decision support system, the newly assimilated simulation states of the real system at time tk are obtained by optimally combining noisy and cluttered measurements with the previously predicted state. Then, the new state estimated is extrapolated to time tk+1 to form the next prediction of the system. The prediction/update cycle repeats as time evolves and new measurements become available. This iterative process is shown in Figure 2, where X represents the abstract state of the system. Accurate assimilation must balance the uncertainty of available measurements with unknown or unpredictable behaviors associated with the model.Figure 2
The relationship between data assimilation and simulation.


## 2.2. DDDS and Data Assimilation

Until now, the simulation models have been using the static information as initial settings to simulate the real systems in the conventional simulation systems. In the dynamic maritime pirate activity modeling, since the static initial parameters fail to capture the real-time environmental changes timely, the simulation states are always quite different from the actual ones of the real system. DDDS has a golden opportunity for development, for the online measurements can reflect the latest updated states of the real system. It means dynamically assimilating online measurements when executing simulations. This idea has been applied in some research areas, such as traffic management, disaster forecasting, environmental science, and wildfire [17, 18]. The simulation system in DDDS paradigm is continually adjusted by assimilating real-time measurements for much more accurate prediction of the studied systems [19]. The rapid development of network and measurement technologies contributes to the higher accuracy and better availability of online measurements, which has facilitated the development of DDDS.DDDS highly depends on data assimilation technology which is in charge of assimilating online measurements into a running simulation system. Data assimilation has the potential to greatly improve the accuracy of simulation systems. As the simulation models become increasingly complicated, if we still make little usage of online measurements from the studied system, traditional simulations will be largely decoupled from real systems. In data assimilation, the measurements are assimilated into the running simulation models to generate more accurate simulation results [20]. The data assimilation algorithms aim to minimize the differences between the simulated states and the real world. There are two main types of data assimilation algorithms: one is the constant statistical algorithms which include 3DVAR [21] and 4DVAR [22] and the other one is the adaptive statistical algorithms, such as Kalman filter [23], ensemble Kalman filter (EnKF) [24], particle filter (PF) [25], and hierarchical Bayesian method [26]. The data assimilation algorithms that have been successfully used in simulation applications include 3DVAR, 4DVAR, Kalman filter, and PF based algorithms.The basic principle of 3DVAR and 4DVAR is to avoid the explicit calculation of the gain matrix and make its inversion by using a minimization procedure of the cost function. They are cost function based and use a suited descent algorithm to iteratively evaluate the cost function. The limitation of variational methods is that they use the minimization of a cost function and they do not provide any estimate of predictive uncertainty. Kalman filter and its variance [27, 28] are important and conventional Bayes estimation techniques used in data assimilation. Kalman filter based data assimilation algorithms are optimal for the assimilation of linear Gaussian systems. The requirements of using Kalman filter based data assimilation algorithms are that the studied systems should be modeled in linear equations and the distributions of processing noises and observation noises are Gaussian. The constraints greatly restrict the applicability of Kalman filter based algorithms in data assimilation for systems that are nonlinear and non-Gaussian. In this paper, we extend the Kalman filter by using RFS based models and unscented transformation.These conventional data assimilation algorithms fail to support data assimilation in the simulation field because of the nonlinear, non-Gaussian properties of a lot of simulation models. The sequential Monte Carlo method is also named as particle filter (PF). PF based data assimilation has been used for some simulation applications including traffic and wildfire. Particle filter is the exact method for the nonlinear and non-Gaussian applications. It has been the most widely used data assimilation algorithm in simulation applications in the recent years [19]. PF based data assimilation in traffic simulation is presented by Xie et al. in [29] and Wu et al. in [30]. PF based data assimilation in wildfire simulation is presented by Xue et al. in [31, 32]. They all use the same nonparametric statistic inference method based on particle filter, because the particle filter based data assimilation algorithms need no assumptions about the distribution or linearity of the studied simulation systems.PF based data assimilation algorithms are widely used, because they are the exact methods for the nonlinear and non-Gaussian systems. But they have the following essential disadvantages:(i)
They are based on the assumption that the studied systems are single dynamic systems that are permanently active. What is more, they use the presumed physical models of individual agents, and the transition of agents from one model to another is less common.(ii)
They are based on the assumption that there is exactly one measurement at each time and the detection is perfect with no false detections or misdetections. In other words, they assume that the sensors never fail to make a measurement of the state and never pick up false measurements.(iii)
They are based on the assumption that the measurements are precise and the measurement models are always precise as well. They always model the measurements as points in the measurement space, and the number and ordering of measurements are previously designated.The standard PF based data assimilation algorithms cannot be used for data assimilation of the dynamic systems that are switching on and off randomly. Switching is quite common in maritime pirate activity, as it did in the case where moving pirate vessels appear and vanish in the designated areas. The restriction in which the studied system must be a single dynamic system limits the applicability of the PF based algorithms. In the maritime pirate activity, multiple dynamic pirate vessels can appear and vanish during the measurement intervals. The measurement models used in PF based data assimilation algorithm are quite unreasonable, since there are always false detections and misdetections in the real world, and the measurement functions are impossible to be precisely known. The PF based algorithm cannot jointly estimate the number of pirates and the states of every pirate either.For the agent based simulation of the maritime pirate activity, we employ RFS to assimilate online measurements into the simulation model. A major task for the data assimilation is to dynamically update the simulation states of the real pirate vessels and then provide the updated states to the simulation model for follow-on simulations. The proposed RFS based data assimilation algorithm can overcome the limitations of the standard PF based data assimilation algorithms effectively.The RFS theory proposed by Goodman et al. as finite set statistics (FISST) is a recently developed theory that unifies much of information fusion under a single Bayesian paradigm [33, 34]. FISST is the first systematic treatment of multisensor multitarget filtering as part of a unified framework for data fusion using random set theory [35, 36]. FISST is a means of directly extending the familiar methodology of single-sensor, single-target Bayes statistics to multisensor multitarget problems. FISST gives some applicable mathematical tools that enable direct application of Bayesian inference. FISST provides the concept of integration/density for set-valued variable. More details can be found in [37].


## 3. RFS Based Models for Data Assimilation

The data assimilation algorithm based on the RFS theory has markedly widened the application fields of data assimilation. It supports both permanently active systems (i.e., present pirates) and on/off systems (i.e., appearing/vanishing pirates). It supports both the single dynamic system and the multiple dynamic systems. It also supports both the perfect detection and the imperfect detection. It can be applied to both precise measurements and the imprecise measurements. It can jointly estimate the number of agents and the state of each of the agents. We will tell how the proposed method overcomes the three limitations of the traditional PF based algorithms under the background of the DDDS of maritime pirate activity.The proposed data assimilation algorithm uses two kinds of system models to describe the real system: one is Markov transition density based simulation model and the other is the measurement likelihood based measurement model. Simulation model is used to describe the state transition of the studied system, and the measurement is used to model the sensors. Both of these models are constructed based on the general data model. The RFS based representation of data (including measurements and states) is the foundation of RFS based data assimilation algorithm and the essential difference between the proposed algorithm and the traditional PF based algorithms.


### 3.1. RFS Based General Data Model

The standard PF based data assimilation assumes precise measurements and measurement models. But the reality is that there is uncertainty in both measurements and measurement models. The uncertainty in the measurements can come from many aspects, including uncertainty due to multiple agents, uncertainty due to unreliability, uncertainty due to false measurements, uncertainty due to deception, uncertainty due to evasive action, and uncertainty due to unknown correlations. The details can be found in [33, 37]. The challenges caused by the uncertainties in the data model have inhibited the maturation of DDDS as a unified, systematic, scientifically founded simulation paradigm. For example, the synthetic aperture radar (SAR) images used in the maritime simulation are statistically uncharacterized and difficult to model. The natural-language statements or rules drawn from knowledge bases are widely used in the maritime pirate activity, and they are unclear and cannot be statistically processed.The limitations of the standard PF based data assimilation algorithms are caused by the fact that they employ vectors to represent the states and measurements. Clearly, the vector based representation cannot represent all occurrences of multiple agent states and, more importantly, does not admit a meaningful and mathematically consistent notion of estimation error. The commonly used vector based assimilation algorithm is summarized in Figure3. We can find that it should assume that the dimension and elements order in each vector are equal. It also needs necessary operations outside of the Bayesian recursion to ensure the consistency of the vectors used for assimilation. The measurement processing step is necessary for vector based algorithms. The determination of newly observed measurement and missed measurement is through vector augmentation or truncation which are very computationally intensive and irreversible. It is also a disadvantage of these algorithms, because it makes the computation complex and unstable. Because of the fundamental inconsistency of the vector based representation, stacking individual measurements into a large vector is not a satisfactory representation.Figure 3
Vector based data assimilation and its mathematical restrictions for a single prediction/update cycle.Here we use RFS to represent the multiple agent states and measurements, because the RFS can represent all possible occurrences of multiple agent states, and distance between sets is a well understood concept. The RFS based data assimilation framework is summarized in Figure4; it enables the joint estimation of errors in cardinality and assimilated states, because it jointly propagates both the estimate of the cardinality of agents and their states. The data association process is automatically dealt with within the FISST framework, and consequently there is no need for data association and measurement management.Figure 4
RFS based data assimilation for a single prediction/update cycle.In the real operational context of maritime, the number of measurements is a random variable, and there is no information on which is the measurement of the state. Actually, the number of measurements is not fixed and the ordering of measurements is not relevant; the measurement is better to be modeled as a RFS. RFS based data model generalizes the state space model for simulation; it takes into account a more realistic situation where the random number of agents and measurements, detection uncertainty, false alarms, and association uncertainty are all taken into consideration.A RFSX=x1,x2,…,xn is a finite set-valued random variable that can be illustrated by a discrete probability distribution ρn=PrX=n and a family of symmetric joint probability densities fnx1,…,xn [33, 38]. The joint probability densities conditioned on cardinality n describe the probability distribution of its elements over the state space. The discrete probability distribution describes the cardinality of RFS X and is defined as follows:(1)ρXn=PrX=n=1n!∫fXy1,…,yndy1⋯dyn.The relationship between the distribution of RFSX and its cardinality is(2)fx1,…,xn=n!·ρn·fnx1,…,xn.We assimilate measurements with a RFS variable formal measurement model that describes sensors’ behaviors and a RFS variable formal simulation model that models the pirates’ behaviors. This means that we use the RFS formulation to treat the collection of individual pirate states as a random finite set-valued state and the collection of individual measurements as a random finite set-valued measurement.


### 3.2. RFS Based Measurement Model

We propose a general detection type measurement model based on the following assumptions as shown in Figure5:(i)
A single sensor with likelihood functiongk+1z∣x=gWk+1z-hk+1x observes a scene involving an unknown number of pirate agents; here hk+1x is the deterministic generation of measurement from state x.(ii)
A single pirate agent with statex generates either a single measurement which means an agent detection, occurring with probability pDx, or a missed detection which means no measurement at all, occurring with probability 1-pDx.(iii)
The false alarm processC is Poisson-distributed in time with expected value λ and distributed in space according to an arbitrary density cz, and ck∣kZk=e-λλk∣kmck∣kz1⋯ck∣kzm is the probability that a set Zk=z1,…,zm of clutter observations will be generated.Figure 5
The major assumptions of RFS based measurement model.For the given predicted multiple agent statesXk=x1,…,xn, the random measurement set collected by the sensor will have the form Zk=ΓkXk∪Ck, where ΓkXk is the agent detection set and it has the form ΓkXk=Γkx1∪⋯∪Γkxn; Γkxj is the detection set for state xj; and Ck is Poisson with mean value λ and spatial distribution cz. We assume that Γkx1,…,Γkxn; Ck is statistically independent.Any given agent can generate either a single measurement or no measurement at all. Consequently,Γkxi must have the form Γkxi=⌀pDxi∩Zi. Here, Zi=hxi,Wi is the sensor-noise model associated with ith state xi; and ⌀p is the discrete random subset of the baseline measurement space O0. Detection uncertainty is modeled as Bernoulli RFS as follows:(3)Pr⌀p=T=1-pD,ifT=⌀,pD·gZi,ifT=O0,0,ifotherwise.The fact that the false alarm processCk is Poisson means that it must have the form Ck=C1,…,CM, where M=Ck is a random nonnegative integer with probability distribution pMm=e-λλm/m!; C1,…,CM are independent, and they are identically distributed random measurement vectors with density cz=gCz.The likelihoodgkZk∣Xk of multiple agents describes the sensor measurement for multiple agents and characterizes the detections, clutters, and agent generated observations. For the state set Xk, the probability density distribution of receiving the measurement set Zk is described by the true likelihood function. The true likelihood function for arbitrary number of agents with missed detections and clutters is (4)gkZk∣Xk=gkΓkXk∪Ck∣Xk=∑W⊆ZkpΓkXkW∣XkpCkZk-W=eλgCZk·gk⌀∣Xk·∑θ∏i:θi>0pDxi·gkzθi∣xi1-pDxi·λczθi;here, gCZk=e-λ∏z∈Zλcz and gk⌀∣Xk=e-λ∏x∈X1-pDx; the summation is taken over all association hypotheses θ:1,…,n→0,1,…,m.


### 3.3. RFS Based Simulation Model

In the agent based simulation of maritime pirate activity, the agents number can be constantly varying because agents can appear and disappear. For example, completely new agents can enter a scene spontaneously, as when a pirate vessel or military vessel suddenly appears in the piracy zones. Existing agents can give rise to new agents through spawning, as when a pirate mothership acts as a floating platform from which smaller pirate vessels are launched and when an apparently single agent at a distance resolves into multiple, closely spaced agents and so on. Agents can likewise leave a scene, as when disappearing from the piracy zones or behind some other occlusion, or they can be damaged or destroyed. All of these instances of multiple agent dynamics should be explicitly modeled by using RFS. If multiple agents are related to birth, death, spawning, merging, and so on, the standard PF based assimilation algorithms are inapplicable, because they fail to accurately estimate the number of agents. In DDDS of maritime pirate activity, we use the unified representation of pirate’s states to get a complete description of the simulation model.We model the time evolution of the multiple agent states (e.g., the states of pirate vessels) by employing RFS. We take the agent birth, death, spawn, and motion into consideration. The RFS based simulation model is based on the following assumptions as shown in Figure6:(i)
The likelihood that a single agent will have state vectorxk at time step k if it had state vector xk-1 at time step k-1 is described by a Markov transition density fk∣k-1xk∣xk-1 corresponding to a single agent simulation model Xk∣k-1=φk-1xk-1,Vk-1. Hereafter this is assumed for notational simplicity to be additive: Xk∣k-1=φk-1xk-1+Vk-1.(ii)
A single agent with statexk-1 at time step k-1 has a probability pSxk=abbr.pSk∣k-1xk-1 of surviving into time step k.(iii)
A single pirate agent with statexk-1 at time step k-1 has a probability density bBk∣k-1Xk-1∣xk-1=abbr.bk∣k-1Bk∣k-1Xk-1∣xk-1 of spawning a set Bk∣k-1Xk-1 of new agents at time step k, as when new pirate vessels launch from a mother ship.(iv)
The probability density that new agents with state setBk will spontaneously appear at time k is bBk=abbr.bk∣k-1Bk, as when some pirate vessels ambush from concealments or enter the zone for the first time.Figure 6
A summary of RFS based simulation model of pirate agents.The multiple agent Markov densityfk∣k-1Xk∣Xk-1 characterizes the time evolution of the multiple agent states and takes the agent births, deaths, and motions into consideration. Here we will give the true multiple agent Markov density function fk∣k-1Xk∣Xk-1 for the proposed simulation model, where Xk=x1,…,xn and Xk-1=x1′,…,xnk-1′. The proposed representation of agents’ states has the mathematical form (5)Xk=Sk∣k-1Xk-1∪Bk∣k-1Xk-1∪Bk,where Xk is the predicted state set, Sk∣k-1Xk-1 are persisting agents that are continuing to chase the merchant, Bk∣k-1Xk-1 are spawned agents that are launched from the mother ship, Bk are spontaneous agents that are ambushing from concealments or enter the designated zone for the first time, Sk∣k-1Xk-1 have the form Sk∣k-1Xk-1=Sk∣k-1x1′∪⋯∪Sk∣k-1xnk-1′, and Sk∣k-1xi′ is the prediction of xi′. Bk∣k-1Xk-1 have the form Bk∣k-1Xk-1=Bk∣k-1x1′∪⋯∪Bk∣k-1xnk-1′; Bk∣k-1xi′ is spawned by xi′ which means that it is launched from the mother ship which has the state of xi′. We also assume that Sk∣k-1x1′,…,Sk∣k-1xnk-1′, Bk∣k-1x1′,…,Bk∣k-1xnk-1′, and Bk are statistically independent.The RFS based simulation model is identical in mathematical form to the RFS based measurement model. The corresponding true multiple agent Markov density is(6)fk∣k-1Xk∣Xk-1=fk∣k-1Sk∣k-1Xk-1∪Bk∣k-1Xk-1∪Bk∣Xk-1=∑W1∩W2=⌀,W2∪W2⊆XkpSk∣k-1W1∣Xk-1pBk∣k-1W2∣Xk-1pBkXk-W1-W2=eμXk-1fBk∣k-1Xk-1Xk·fk∣k-1⌀∣Xk-1·∑θ∏i:θi>0pSxi′·fk∣k-1xθi∣xi′1-pSxi′·μXk-1bxθi∣Xk-1.The summation is taken over all association hypothesesθ:1,…,nk-1→0,1,…,nk, where (7)fk∣k-1⌀∣Xk-1=e-μXk-1∏xk-1∈Xk-11-pSxk-1,μXk-1=μ0+μ1x′+⋯+μnk-1x′,bx∣Xk-1=μ0b0x+μ1x′bx∣x1′+⋯+μnk-1x′bx∣xnk-1′μ0+μ1x′+⋯+μnk-1x′.Here,μ0 is the expected number of spontaneously generated new agents and b0x is their physical distribution. Also, μxi′ is the expected number of new agents spawned by an agent with previous state xi′, and bx∣xi′ is their physical distribution.The calculation of the true Markov density is impractical. Here, we use another new variable to deal with the complex computation. The difficulty of RFS based Bayesian inference for multiple agents is its computational complexity. To solve this problem, Mahler proposed the PHD (probability hypothesis density) filter [34]. The PHD Dkkx∣Zk of the posterior probability density fkkXk∣Zk is a probability density function defined on the single agent state x∈X0 as follows:(8)Dkk=abbr.Dk∣kx∣Zk≜∫fk∣kx∪Xk∣ZkδXk.In point theory,Dk∣kx is defined as the intensity density [39]. It is not a probability density. Dk∣kx represents the density of expected number of points at x. Given any region S of single pirate agent state space X0, the integral ∫SDk∣kxdx=vk∣k is the expected number of pirate agents in S, where vk∣k∈R is the expected number of objects in the state space S. In particular, if S=X0 is the entire state space, then Nk∣k≜∫Dk∣kxdx is the total expected number of pirate agents in the scene.Here, we assume that then PHD of the RFSBk for the new born pirate agents is a Gaussian mixture defined as follows: (9)γkx=∑i=1Jγ,kwγ,kiNx;mγ,ki,Pγ,ki,where wγ,ki, mγ,ki, and Pγ,ki are the weights, means, and covariance matrices of the mixture birth PHD.


## 3.1. RFS Based General Data Model

The standard PF based data assimilation assumes precise measurements and measurement models. But the reality is that there is uncertainty in both measurements and measurement models. The uncertainty in the measurements can come from many aspects, including uncertainty due to multiple agents, uncertainty due to unreliability, uncertainty due to false measurements, uncertainty due to deception, uncertainty due to evasive action, and uncertainty due to unknown correlations. The details can be found in [33, 37]. The challenges caused by the uncertainties in the data model have inhibited the maturation of DDDS as a unified, systematic, scientifically founded simulation paradigm. For example, the synthetic aperture radar (SAR) images used in the maritime simulation are statistically uncharacterized and difficult to model. The natural-language statements or rules drawn from knowledge bases are widely used in the maritime pirate activity, and they are unclear and cannot be statistically processed.The limitations of the standard PF based data assimilation algorithms are caused by the fact that they employ vectors to represent the states and measurements. Clearly, the vector based representation cannot represent all occurrences of multiple agent states and, more importantly, does not admit a meaningful and mathematically consistent notion of estimation error. The commonly used vector based assimilation algorithm is summarized in Figure3. We can find that it should assume that the dimension and elements order in each vector are equal. It also needs necessary operations outside of the Bayesian recursion to ensure the consistency of the vectors used for assimilation. The measurement processing step is necessary for vector based algorithms. The determination of newly observed measurement and missed measurement is through vector augmentation or truncation which are very computationally intensive and irreversible. It is also a disadvantage of these algorithms, because it makes the computation complex and unstable. Because of the fundamental inconsistency of the vector based representation, stacking individual measurements into a large vector is not a satisfactory representation.Figure 3
Vector based data assimilation and its mathematical restrictions for a single prediction/update cycle.Here we use RFS to represent the multiple agent states and measurements, because the RFS can represent all possible occurrences of multiple agent states, and distance between sets is a well understood concept. The RFS based data assimilation framework is summarized in Figure4; it enables the joint estimation of errors in cardinality and assimilated states, because it jointly propagates both the estimate of the cardinality of agents and their states. The data association process is automatically dealt with within the FISST framework, and consequently there is no need for data association and measurement management.Figure 4
RFS based data assimilation for a single prediction/update cycle.In the real operational context of maritime, the number of measurements is a random variable, and there is no information on which is the measurement of the state. Actually, the number of measurements is not fixed and the ordering of measurements is not relevant; the measurement is better to be modeled as a RFS. RFS based data model generalizes the state space model for simulation; it takes into account a more realistic situation where the random number of agents and measurements, detection uncertainty, false alarms, and association uncertainty are all taken into consideration.A RFSX=x1,x2,…,xn is a finite set-valued random variable that can be illustrated by a discrete probability distribution ρn=PrX=n and a family of symmetric joint probability densities fnx1,…,xn [33, 38]. The joint probability densities conditioned on cardinality n describe the probability distribution of its elements over the state space. The discrete probability distribution describes the cardinality of RFS X and is defined as follows:(1)ρXn=PrX=n=1n!∫fXy1,…,yndy1⋯dyn.The relationship between the distribution of RFSX and its cardinality is(2)fx1,…,xn=n!·ρn·fnx1,…,xn.We assimilate measurements with a RFS variable formal measurement model that describes sensors’ behaviors and a RFS variable formal simulation model that models the pirates’ behaviors. This means that we use the RFS formulation to treat the collection of individual pirate states as a random finite set-valued state and the collection of individual measurements as a random finite set-valued measurement.


## 3.2. RFS Based Measurement Model

We propose a general detection type measurement model based on the following assumptions as shown in Figure5:(i)
A single sensor with likelihood functiongk+1z∣x=gWk+1z-hk+1x observes a scene involving an unknown number of pirate agents; here hk+1x is the deterministic generation of measurement from state x.(ii)
A single pirate agent with statex generates either a single measurement which means an agent detection, occurring with probability pDx, or a missed detection which means no measurement at all, occurring with probability 1-pDx.(iii)
The false alarm processC is Poisson-distributed in time with expected value λ and distributed in space according to an arbitrary density cz, and ck∣kZk=e-λλk∣kmck∣kz1⋯ck∣kzm is the probability that a set Zk=z1,…,zm of clutter observations will be generated.Figure 5
The major assumptions of RFS based measurement model.For the given predicted multiple agent statesXk=x1,…,xn, the random measurement set collected by the sensor will have the form Zk=ΓkXk∪Ck, where ΓkXk is the agent detection set and it has the form ΓkXk=Γkx1∪⋯∪Γkxn; Γkxj is the detection set for state xj; and Ck is Poisson with mean value λ and spatial distribution cz. We assume that Γkx1,…,Γkxn; Ck is statistically independent.Any given agent can generate either a single measurement or no measurement at all. Consequently,Γkxi must have the form Γkxi=⌀pDxi∩Zi. Here, Zi=hxi,Wi is the sensor-noise model associated with ith state xi; and ⌀p is the discrete random subset of the baseline measurement space O0. Detection uncertainty is modeled as Bernoulli RFS as follows:(3)Pr⌀p=T=1-pD,ifT=⌀,pD·gZi,ifT=O0,0,ifotherwise.The fact that the false alarm processCk is Poisson means that it must have the form Ck=C1,…,CM, where M=Ck is a random nonnegative integer with probability distribution pMm=e-λλm/m!; C1,…,CM are independent, and they are identically distributed random measurement vectors with density cz=gCz.The likelihoodgkZk∣Xk of multiple agents describes the sensor measurement for multiple agents and characterizes the detections, clutters, and agent generated observations. For the state set Xk, the probability density distribution of receiving the measurement set Zk is described by the true likelihood function. The true likelihood function for arbitrary number of agents with missed detections and clutters is (4)gkZk∣Xk=gkΓkXk∪Ck∣Xk=∑W⊆ZkpΓkXkW∣XkpCkZk-W=eλgCZk·gk⌀∣Xk·∑θ∏i:θi>0pDxi·gkzθi∣xi1-pDxi·λczθi;here, gCZk=e-λ∏z∈Zλcz and gk⌀∣Xk=e-λ∏x∈X1-pDx; the summation is taken over all association hypotheses θ:1,…,n→0,1,…,m.


## 3.3. RFS Based Simulation Model

In the agent based simulation of maritime pirate activity, the agents number can be constantly varying because agents can appear and disappear. For example, completely new agents can enter a scene spontaneously, as when a pirate vessel or military vessel suddenly appears in the piracy zones. Existing agents can give rise to new agents through spawning, as when a pirate mothership acts as a floating platform from which smaller pirate vessels are launched and when an apparently single agent at a distance resolves into multiple, closely spaced agents and so on. Agents can likewise leave a scene, as when disappearing from the piracy zones or behind some other occlusion, or they can be damaged or destroyed. All of these instances of multiple agent dynamics should be explicitly modeled by using RFS. If multiple agents are related to birth, death, spawning, merging, and so on, the standard PF based assimilation algorithms are inapplicable, because they fail to accurately estimate the number of agents. In DDDS of maritime pirate activity, we use the unified representation of pirate’s states to get a complete description of the simulation model.We model the time evolution of the multiple agent states (e.g., the states of pirate vessels) by employing RFS. We take the agent birth, death, spawn, and motion into consideration. The RFS based simulation model is based on the following assumptions as shown in Figure6:(i)
The likelihood that a single agent will have state vectorxk at time step k if it had state vector xk-1 at time step k-1 is described by a Markov transition density fk∣k-1xk∣xk-1 corresponding to a single agent simulation model Xk∣k-1=φk-1xk-1,Vk-1. Hereafter this is assumed for notational simplicity to be additive: Xk∣k-1=φk-1xk-1+Vk-1.(ii)
A single agent with statexk-1 at time step k-1 has a probability pSxk=abbr.pSk∣k-1xk-1 of surviving into time step k.(iii)
A single pirate agent with statexk-1 at time step k-1 has a probability density bBk∣k-1Xk-1∣xk-1=abbr.bk∣k-1Bk∣k-1Xk-1∣xk-1 of spawning a set Bk∣k-1Xk-1 of new agents at time step k, as when new pirate vessels launch from a mother ship.(iv)
The probability density that new agents with state setBk will spontaneously appear at time k is bBk=abbr.bk∣k-1Bk, as when some pirate vessels ambush from concealments or enter the zone for the first time.Figure 6
A summary of RFS based simulation model of pirate agents.The multiple agent Markov densityfk∣k-1Xk∣Xk-1 characterizes the time evolution of the multiple agent states and takes the agent births, deaths, and motions into consideration. Here we will give the true multiple agent Markov density function fk∣k-1Xk∣Xk-1 for the proposed simulation model, where Xk=x1,…,xn and Xk-1=x1′,…,xnk-1′. The proposed representation of agents’ states has the mathematical form (5)Xk=Sk∣k-1Xk-1∪Bk∣k-1Xk-1∪Bk,where Xk is the predicted state set, Sk∣k-1Xk-1 are persisting agents that are continuing to chase the merchant, Bk∣k-1Xk-1 are spawned agents that are launched from the mother ship, Bk are spontaneous agents that are ambushing from concealments or enter the designated zone for the first time, Sk∣k-1Xk-1 have the form Sk∣k-1Xk-1=Sk∣k-1x1′∪⋯∪Sk∣k-1xnk-1′, and Sk∣k-1xi′ is the prediction of xi′. Bk∣k-1Xk-1 have the form Bk∣k-1Xk-1=Bk∣k-1x1′∪⋯∪Bk∣k-1xnk-1′; Bk∣k-1xi′ is spawned by xi′ which means that it is launched from the mother ship which has the state of xi′. We also assume that Sk∣k-1x1′,…,Sk∣k-1xnk-1′, Bk∣k-1x1′,…,Bk∣k-1xnk-1′, and Bk are statistically independent.The RFS based simulation model is identical in mathematical form to the RFS based measurement model. The corresponding true multiple agent Markov density is(6)fk∣k-1Xk∣Xk-1=fk∣k-1Sk∣k-1Xk-1∪Bk∣k-1Xk-1∪Bk∣Xk-1=∑W1∩W2=⌀,W2∪W2⊆XkpSk∣k-1W1∣Xk-1pBk∣k-1W2∣Xk-1pBkXk-W1-W2=eμXk-1fBk∣k-1Xk-1Xk·fk∣k-1⌀∣Xk-1·∑θ∏i:θi>0pSxi′·fk∣k-1xθi∣xi′1-pSxi′·μXk-1bxθi∣Xk-1.The summation is taken over all association hypothesesθ:1,…,nk-1→0,1,…,nk, where (7)fk∣k-1⌀∣Xk-1=e-μXk-1∏xk-1∈Xk-11-pSxk-1,μXk-1=μ0+μ1x′+⋯+μnk-1x′,bx∣Xk-1=μ0b0x+μ1x′bx∣x1′+⋯+μnk-1x′bx∣xnk-1′μ0+μ1x′+⋯+μnk-1x′.Here,μ0 is the expected number of spontaneously generated new agents and b0x is their physical distribution. Also, μxi′ is the expected number of new agents spawned by an agent with previous state xi′, and bx∣xi′ is their physical distribution.The calculation of the true Markov density is impractical. Here, we use another new variable to deal with the complex computation. The difficulty of RFS based Bayesian inference for multiple agents is its computational complexity. To solve this problem, Mahler proposed the PHD (probability hypothesis density) filter [34]. The PHD Dkkx∣Zk of the posterior probability density fkkXk∣Zk is a probability density function defined on the single agent state x∈X0 as follows:(8)Dkk=abbr.Dk∣kx∣Zk≜∫fk∣kx∪Xk∣ZkδXk.In point theory,Dk∣kx is defined as the intensity density [39]. It is not a probability density. Dk∣kx represents the density of expected number of points at x. Given any region S of single pirate agent state space X0, the integral ∫SDk∣kxdx=vk∣k is the expected number of pirate agents in S, where vk∣k∈R is the expected number of objects in the state space S. In particular, if S=X0 is the entire state space, then Nk∣k≜∫Dk∣kxdx is the total expected number of pirate agents in the scene.Here, we assume that then PHD of the RFSBk for the new born pirate agents is a Gaussian mixture defined as follows: (9)γkx=∑i=1Jγ,kwγ,kiNx;mγ,ki,Pγ,ki,where wγ,ki, mγ,ki, and Pγ,ki are the weights, means, and covariance matrices of the mixture birth PHD.


## 4. Agent Model of Vessels

Multiple agent based simulation models are more suitable for characterizing the complex systems. They can capture the properties of reasoning processes. The agent based modeling methods are perfect for the simulation of maritime pirate activity. Because there are a lot of priorities, objectives, constraints, and complex scheduling dependencies in the actual maritime pirate activities, the agent based simulation has been applied to solve these problems successfully.Here, the agent based simulation models are used to represent the behaviors of the real pirate vessels. The pirate vessels are characterized by autonomous agents. These agent models describe the movement and interactions of different kinds of pirate vessels. If there are no interactions between the vessel agents, each vessel agent pursues its individual goals. But just like the real pirate vessels, the pirate vessel agents can also interact with other kinds of vessels. And these interactions are difficult to model by using the mathematical models or DEVS.The data assimilation algorithm is implemented in the merchant vessel agent model. By assimilating the measurements from the real system, the merchant vessel agent can interact with the real system more conveniently and effectively. The interactions between the merchant agent and the real system are shown in Figure7.Figure 7
Interactions between merchant vessel agent and real system.


## 5. RFS Based Data Assimilation Algorithm

We propose using a new RFS based data assimilation algorithm by combining the basic ideas of Gaussian mixture approximation, Cardinalized Probability Hypothesis Density (CPHD) filter, and the unscented Kalman filter [40] together. Unlike the standard unscented Kalman filter, the proposed algorithm represents the states by RFS rather than by vector.LetpkX∣Z1:k denote the multiple agent posterior probability density. Since we have got the RFS based measurement model and simulation model, we can use these two models to get the optimal Bayes filter. The optimal Bayesian inference is captured by the recursion as follows:(10)pk∣k-1Xk∣Z1:k-1=∫fk∣k-1Xk∣Xpk-1X∣Z1:k-1μSdX,pkXk∣Z1:k=gkZk∣Xkpk∣k-1Xk∣Z1:k-1∫gkZk∣Xpk∣k-1X∣Z1:k-1μSdX,where μS is the reference measure on all finite subsets of the measurement space [41].The CPHD filter is introduced by Mahler in [42], and it propagates both the PHD and the discrete distribution of the pirates’ number. The discrete distribution of pirate numbers can be arbitrary.LetDk∣kx∣Z1:k=abbr.Dk∣kx denote the PHD associated with the posterior multiple agent states, and let ρk∣kn∣Z1:n=abbr.ρk∣kn denote the posterior distribution of pirates’ number at time k. Also, let ρΓ,k∣k-1n=abbr.ρΓn denote the distribution of new agents’ number at time k. The proposed data assimilation algorithm propagates both the cardinality distribution ρk∣kn and the PHD Dk∣kx.The main idea of the proposed data assimilation algorithm is shown in Figure8. It is a rollback based algorithm. In this algorithm, the Gaussian components are recursively generated by the rollback based sampling step. The rollback operation means reversing or undoing something. The rollback in this paper has its natural meaning which means restoring the simulation to a previous state. The agent based simulation system is used in the prediction step of the data assimilation algorithm. The rollback strategy is used in the simulation running during the prediction step. In the prediction step, the agent based simulation system is the parallel system of the real system. The update step is the central step of the data assimilation algorithm. The update step uses the predicted measurement statistics, real-time measurements, and predicted cardinality as inputs.Figure 8
Gaussian mixture based implementation of RFS based data assimilation algorithm.


### 5.1. Predict with RFS Based Simulation Model

If the number of pirate agents that continue to attack the merchant agent isj, the predicted cardinality of these survived agents will be(11)ρS,k∣k-1j=∑l=j∞CjlpS,Dk-1k-1j1-pS,Dk-1∣k-1l-j1,Dk-1∣k-1ρk-1∣k-1l=∑l=j∞l!j!l-j!pSj1-pSl-jρk-1∣k-1l.Then the predicted distribution of pirates’ number can be characterized by the following convolution:(12)ρk∣k-1n=∑j=0nρΓ,kn-jρS,k∣k-1Dk-1∣k-1,ρk-1∣k-1j=∑j=0nρS,k∣k-1jρΓn-j.Based on the RFS based simulation model of pirate agent movement, the predicted PHD is given as follows:(13)Dk∣k-1x=γk∣k-1x+pS∫fk∣k-1x∣x′Dk-1∣k-1x′dx′,where γk∣k-1x denotes the PHD of the RFS of new agents appearing between time k-1 and timek.Since the Markov transition density of the single vessel agent is nonlinear, here we use the unscented transformation to approximate the moments offk∣k-1x∣x′ since we know it can be transformed into a linear-Gaussian model by using the unscented transformation.Under the assumption that the posterior PHD at timek-1 is the Gaussian mixture denoted by Dk-1∣k-1x=∑i=1Jk-1wk-1iNx;mk-1i,Pk-1i, the predicted PHD Dk∣k-1x is also a Gaussian mixture under the assumption Dk∣k-1x=DS,k∣k-1x+γkx, where γkx is given as in (9) and the predicted PHD of survival agents is(14)DS,k∣k-1x=pS∫fk∣k-1x∣x′Dk-1∣k-1x′dx′=pS∑j=1Jk-1wk-1jNx;mS,k∣k-1j,PS,k∣k-1j.For each Gaussian component, we should compute the predicted mean valuemS,k∣k-1j and covariance PS,k∣k-1j at time k by using the unscented transformation and the simulation model. For the ith Gaussian component, its predicted mean value mS,k∣k-1j and predicted covariance PS,k∣k-1j are computed by the following steps. Here, we assume that the mean of the ith Gaussian component is μk-1j and its covariance is Ck-1i.Step 1.
Choose the set of sigma points and their corresponding weights denoted byxk-1l,ull=0L to approximate mean μki and covariance Cki as introduced in [43].Step 2.
Run the agent based simulation system to simulation timetk and get the transformed sigma points denoted by xk∣k-1l:=φkxk-1l,vk-1l,l=0,…,L.Step 3.
Compute the predicted state statistics:(15)mk∣k-1i=∑l=0Lulxk∣k-1l,Pk∣k-1i=∑l=0Lulxk∣k-1l-mk∣k-1ixk∣k-1l-mk∣k-1iT.Step 4.
Choose the sigma pointsxk∣k-11,…,xk∣k-1L and their corresponding weights u1,…,uL to approximate mean mk∣k-1i and covariance Pk∣k-1i.Step 5.
Generate the transformed sigma pointszk∣k-1l:=hkxk∣k-1l,εkl,l=0,…,L.Step 6.
Calculate the predicted measurements and their corresponding characteristic:(16)ηk∣k-1i=∑l=0Lulzk∣k-1l,Ski=∑l=0Lulzk∣k-1l-ηk∣k-1izk∣k-1l-ηk∣k-1iT,Gki=∑l=0Lulxk∣k-1l-mk∣k-1izk∣k-1l-ηk∣k-1iT.Step 7.
Compute the Kalman gain and posterior covariance matrix:(17)Kki=GkiSki-1,Pk∣ki=Pk∣k-1i-GkiSki-1GkiT.


### 5.2. Update with Online Measurements

Without measurements, the covariance matrix grows exponentially over time. Assimilating the real-time measurements can be regarded as calibrating or estimating the real-time state of the simulation. The update step of the proposed data assimilation algorithm is depending on the Bayesian inference. The updated cardinality is as follows:(18)ρk∣kn=Υk,Dk∣k-1;Zk0nρk∣k-1nΥk,Dk∣k-1;Zk0,ρk∣k-1.Under the assumption that the predicted PHDDk∣k-1x at time k is the Gaussian mixture denoted by Dk∣k-1x=∑i=1Jk∣k-1wk∣k-1iNx;mk∣k-1i,Pk∣k-1i, the updated PHD Dk∣kx can also be represented by the Gaussian mixture as follows:(19)Dk∣kx=gkZk∣XkDk∣k-1x∫gkZk∣XkDk∣k-1xdx≅∑z∈ZkpDgkz∣x∫pDgkz∣xDk∣k-1xdx+cz+1-pDDk∣k-1x=∑z∈ZkΥk,Dk∣k-1;Zk∖z1,ρk∣k-1Υk,Dk∣k-1;Zk0,ρk∣k-1pD1,ckgkz∣xckz+1-pDΥk,Dk∣k-1;Zk1,ρk∣k-1Υk,Dk∣k-1;Zk0,ρk∣k-1Dk∣k-1x=1-pDΥk,Dk∣k-1;Zk1,ρk∣k-1Υk,Dk∣k-1;Zk0,ρk∣k-1Dk∣k-1x+∑z∈Zk∑j=1Jk∣k-1wkjzNx;mkjz,Pkj,where (20)wk∣k-1=wk∣k-11,…,wk∣k-1Jk∣k-1T,qkz=qk1z,…,qkJk∣k-1zT,qkjz=Nz;ηk∣k-1i,Sk∣k-1i,wklJk∣k-1+j=pDwk∣k-1jqkjzΥk,Dk∣k-1;Zk∖z1,ρk∣k-1Υk,Dk∣k-1;Zk0,ρk∣k-1gkz∣xcz,mkj=mk∣k-1j+Kkjz-ηk∣k-1j.ρ
1
,
ρ
2 denotes the inner product defined by ρ1,ρ2=∑n=0∞ρ1nρ2n;gk·∣x is the single pirate agent measurement likelihood density at time k based on state x; ck· is the PHD of clutter measurements at time k. The sequence Υk,D,Znun is defined for u∈0,1 as follows:(21)Υk,D,Znun=∑j=0min⁡Z,nZ-j!ρc,kZ-j·Pj+un1-pD∫Dxdx-j+u·ejΞkD,Z.The termρc,kn in (21) is the distribution of the numbers of clutters at time k. The other terms are given as follows:(22)ΞkD,Z=pD∫Dxgkz∣xdxcz:z∈Z,ejZ=∑W⊆Z,W=j∏ξ∈Wξ,Pln=n!n-l!,where ejZ is the elementary symmetric function of order j for the random finite set Z; the details can be found in [37].Letn^ be the integer nearest in value to ρk∣kn. Then the n^ highest peaks of Dk∣ky,0 are corresponding to the states of the agents. That is, let y^1,…,y^n^ be such that Dk∣ky^1,0,…,Dk∣ky^n^,0 are the n^ largest values of Dk∣k·,0. Then Y^k∣k=y^1,…,y^n^ is an estimate of the multiple agent states. The assimilated agent number is n^k∣k=arg⁡sup⁡ρk∣kn∣Zk by using MAP estimation.


## 5.1. Predict with RFS Based Simulation Model

If the number of pirate agents that continue to attack the merchant agent isj, the predicted cardinality of these survived agents will be(11)ρS,k∣k-1j=∑l=j∞CjlpS,Dk-1k-1j1-pS,Dk-1∣k-1l-j1,Dk-1∣k-1ρk-1∣k-1l=∑l=j∞l!j!l-j!pSj1-pSl-jρk-1∣k-1l.Then the predicted distribution of pirates’ number can be characterized by the following convolution:(12)ρk∣k-1n=∑j=0nρΓ,kn-jρS,k∣k-1Dk-1∣k-1,ρk-1∣k-1j=∑j=0nρS,k∣k-1jρΓn-j.Based on the RFS based simulation model of pirate agent movement, the predicted PHD is given as follows:(13)Dk∣k-1x=γk∣k-1x+pS∫fk∣k-1x∣x′Dk-1∣k-1x′dx′,where γk∣k-1x denotes the PHD of the RFS of new agents appearing between time k-1 and timek.Since the Markov transition density of the single vessel agent is nonlinear, here we use the unscented transformation to approximate the moments offk∣k-1x∣x′ since we know it can be transformed into a linear-Gaussian model by using the unscented transformation.Under the assumption that the posterior PHD at timek-1 is the Gaussian mixture denoted by Dk-1∣k-1x=∑i=1Jk-1wk-1iNx;mk-1i,Pk-1i, the predicted PHD Dk∣k-1x is also a Gaussian mixture under the assumption Dk∣k-1x=DS,k∣k-1x+γkx, where γkx is given as in (9) and the predicted PHD of survival agents is(14)DS,k∣k-1x=pS∫fk∣k-1x∣x′Dk-1∣k-1x′dx′=pS∑j=1Jk-1wk-1jNx;mS,k∣k-1j,PS,k∣k-1j.For each Gaussian component, we should compute the predicted mean valuemS,k∣k-1j and covariance PS,k∣k-1j at time k by using the unscented transformation and the simulation model. For the ith Gaussian component, its predicted mean value mS,k∣k-1j and predicted covariance PS,k∣k-1j are computed by the following steps. Here, we assume that the mean of the ith Gaussian component is μk-1j and its covariance is Ck-1i.Step 1.
Choose the set of sigma points and their corresponding weights denoted byxk-1l,ull=0L to approximate mean μki and covariance Cki as introduced in [43].Step 2.
Run the agent based simulation system to simulation timetk and get the transformed sigma points denoted by xk∣k-1l:=φkxk-1l,vk-1l,l=0,…,L.Step 3.
Compute the predicted state statistics:(15)mk∣k-1i=∑l=0Lulxk∣k-1l,Pk∣k-1i=∑l=0Lulxk∣k-1l-mk∣k-1ixk∣k-1l-mk∣k-1iT.Step 4.
Choose the sigma pointsxk∣k-11,…,xk∣k-1L and their corresponding weights u1,…,uL to approximate mean mk∣k-1i and covariance Pk∣k-1i.Step 5.
Generate the transformed sigma pointszk∣k-1l:=hkxk∣k-1l,εkl,l=0,…,L.Step 6.
Calculate the predicted measurements and their corresponding characteristic:(16)ηk∣k-1i=∑l=0Lulzk∣k-1l,Ski=∑l=0Lulzk∣k-1l-ηk∣k-1izk∣k-1l-ηk∣k-1iT,Gki=∑l=0Lulxk∣k-1l-mk∣k-1izk∣k-1l-ηk∣k-1iT.Step 7.
Compute the Kalman gain and posterior covariance matrix:(17)Kki=GkiSki-1,Pk∣ki=Pk∣k-1i-GkiSki-1GkiT.


## 5.2. Update with Online Measurements

Without measurements, the covariance matrix grows exponentially over time. Assimilating the real-time measurements can be regarded as calibrating or estimating the real-time state of the simulation. The update step of the proposed data assimilation algorithm is depending on the Bayesian inference. The updated cardinality is as follows:(18)ρk∣kn=Υk,Dk∣k-1;Zk0nρk∣k-1nΥk,Dk∣k-1;Zk0,ρk∣k-1.Under the assumption that the predicted PHDDk∣k-1x at time k is the Gaussian mixture denoted by Dk∣k-1x=∑i=1Jk∣k-1wk∣k-1iNx;mk∣k-1i,Pk∣k-1i, the updated PHD Dk∣kx can also be represented by the Gaussian mixture as follows:(19)Dk∣kx=gkZk∣XkDk∣k-1x∫gkZk∣XkDk∣k-1xdx≅∑z∈ZkpDgkz∣x∫pDgkz∣xDk∣k-1xdx+cz+1-pDDk∣k-1x=∑z∈ZkΥk,Dk∣k-1;Zk∖z1,ρk∣k-1Υk,Dk∣k-1;Zk0,ρk∣k-1pD1,ckgkz∣xckz+1-pDΥk,Dk∣k-1;Zk1,ρk∣k-1Υk,Dk∣k-1;Zk0,ρk∣k-1Dk∣k-1x=1-pDΥk,Dk∣k-1;Zk1,ρk∣k-1Υk,Dk∣k-1;Zk0,ρk∣k-1Dk∣k-1x+∑z∈Zk∑j=1Jk∣k-1wkjzNx;mkjz,Pkj,where (20)wk∣k-1=wk∣k-11,…,wk∣k-1Jk∣k-1T,qkz=qk1z,…,qkJk∣k-1zT,qkjz=Nz;ηk∣k-1i,Sk∣k-1i,wklJk∣k-1+j=pDwk∣k-1jqkjzΥk,Dk∣k-1;Zk∖z1,ρk∣k-1Υk,Dk∣k-1;Zk0,ρk∣k-1gkz∣xcz,mkj=mk∣k-1j+Kkjz-ηk∣k-1j.ρ
1
,
ρ
2 denotes the inner product defined by ρ1,ρ2=∑n=0∞ρ1nρ2n;gk·∣x is the single pirate agent measurement likelihood density at time k based on state x; ck· is the PHD of clutter measurements at time k. The sequence Υk,D,Znun is defined for u∈0,1 as follows:(21)Υk,D,Znun=∑j=0min⁡Z,nZ-j!ρc,kZ-j·Pj+un1-pD∫Dxdx-j+u·ejΞkD,Z.The termρc,kn in (21) is the distribution of the numbers of clutters at time k. The other terms are given as follows:(22)ΞkD,Z=pD∫Dxgkz∣xdxcz:z∈Z,ejZ=∑W⊆Z,W=j∏ξ∈Wξ,Pln=n!n-l!,where ejZ is the elementary symmetric function of order j for the random finite set Z; the details can be found in [37].Letn^ be the integer nearest in value to ρk∣kn. Then the n^ highest peaks of Dk∣ky,0 are corresponding to the states of the agents. That is, let y^1,…,y^n^ be such that Dk∣ky^1,0,…,Dk∣ky^n^,0 are the n^ largest values of Dk∣k·,0. Then Y^k∣k=y^1,…,y^n^ is an estimate of the multiple agent states. The assimilated agent number is n^k∣k=arg⁡sup⁡ρk∣kn∣Zk by using MAP estimation.


## 6. Experiments

### 6.1. Experimental Design

We take advantage of the identical-twin experiment introduced in [31, 32] to verify the proposed data assimilation algorithm in ideal situations. In the identical-twin experiment, the simulation whose results are considered as real states is first executed with corresponding measurements which are collected during the process. And the measurements here are treated as the real measurements, since they are collected from the real states. Accordingly, we assimilate the measurements into DDDS by using the proposed data assimilation algorithm and then check whether the assimilated results are close to the real system states or not.We utilize three terms in the experiments, real pirate movement, assimilated pirate movement, and simulated pirate movement, to represent the results of the experiments. A real pirate movement is the simulated pirate movement from which the online measurements are obtained. A simulated pirate movement is the simulation states based on some erroneous parameters, for example, imprecise turn rate, which aims to comply with the fact that pirate movement simulation usually relies on imprecise parameters compared to real pirate movement. Finally, an assimilated pirate movement is the data assimilation enhanced simulation states based on the same erroneous parameters as in the simulated pirate movement. The goal of these experiments is to verify that an assimilated pirate movement creates more accurate simulation states by assimilating online measurements from the real pirate movement, even if it uses the erroneous parameters in the simulated pirate movement.


### 6.2. Experiments on Single Pirate Vessel in Clutter

The state space for the pirate vessel is four-dimensional with two position dimensions in the normalx,y Cartesian coordinate system and two velocity coordinates, x˙,y˙. We establish the 2D horizontal turn motion model relying on ship kinematics. We use a widely used ship motion model [44]. The dynamics equation is as follows:(23)xk=xk-1+Tvk-1sin⁡ϕk-1,yk=yk-1+Tvk-1cos⁡ϕk-1,ϕk=ϕk-1+Tvk-1Ωk-1+12Ωk-1-Ω0Tvk-1τeTvk-1τ,Ωk=Ωk-1eTvk-1τ+Ω01-eTvk-1τ,vk-1=K·v0=v01+1.9Ωk-12L2-1;here, τ=-p/2+p2/4-q/L, Ω0=Ω/v0, and x,y, Ω, and ϕ are ship position, velocity vector turn rate, and heading, respectively. T is the sampling interval.The simulation modelxk=φkxk-1,vk-1 of the vessel motion is nonlinear, where φk is the nonlinear state transition functions and vk-1 is an independent zero-mean Gaussian process noise; and Qk-1 is the covariance matrix. The vessel simulation model can be reformulated as(24)xk=φk-1xk-1+G·ϕ˙k-1,ϕ˙k=ϕ˙k-1+T·uk-1,where G=T2/2T0000T2/2TT,ϕ˙k-1~N·;0,σϕ˙2I,uk-1~N·;0;σu2I, T = 1 s, σϕ˙ = 15 m/s2, and σu=π/180 rad/s. The probability of survival is pS,kx = 0.99. The probability of detection is pD,kx = 0.98. The speed of the pirate vessel is 25–30 knots.The sensor measurement model in the experiments is as follows:(25)zk=arctan⁡xkykxk2+yk2+εk,where measurement noise εk~N·;0,Rk with Rk=diag(σθ2,σr2T); in the simulation experiments, we assume that σθ = 2(π/180) rad, and σr = 10 m.The false detection obeys the Poisson RFS with PHDλc = 3.2 × 10−3 (rad·m)−1 over the disc [0, 2π] rad × [0,2000] m.In our experiments, the differences between a real pirate movement and a simulated one are caused by the imprecise turn rateϕ˙. We use the imprecise turn rate as the erroneous parameters. The real turn rate is π/60≈0.0524 rad/s with random variances added every second. Our first experiment case uses the real turn rate to run the agent based simulation system and collects the measurements at every time step. Our second experiment case introduces errors to the turn rate and makes the other configurations exactly the same as the first experiment. The erroneous turn rate is 0.0524+0.01/3≈0.05573 rad/s with the same random variances in the first experiment. Our third experiment case uses the erroneous turn rate as the second experiment and the other configurations remain the same as the other two experiments. In this experiment case, the measurements collected in the first experiment are used as well as the proposed data assimilation algorithm.We assume that the merchant vessel is stationary at a particular location, which is the center in the Cartesian coordinates. In this way, we can get the pirate’s states. The states of the pirate vessel generated by the three experiments are shown in Figure9.Figure 9
States of identical-twin experiments on single pirate vessel in clutter.This set of experiments verifies the effectiveness of data assimilation algorithm when the turn rate has errors. From Figure9, we know that the real pirate vessel movement and the simulated one have large deviations due to the turn rate errors. By assimilating the measurements, the assimilated pirate vessel movement is obtained. From this figure, it is observed that the assimilated pirate vessel movement is much closer to the real one than the simulated one. Specifically, due to the erroneous turn rate, the simulated pirate states gradually deviate from the real pirate states. Using the proposed data assimilation algorithm, the assimilated pirate states overcome this problem and match the real pirate states with minor differences.Figure10(a) shows the PHD of positions and Figure 10(b) shows the contours of the marginal PHD in all time steps. We can learn that the PHD is gradually improved as the simulation time evolves. We can find that, in the initial phase, the uncertainty of the PHD is large due to the limited and imprecise prior information, and as the simulation time evolves, the uncertainty gradually decreases. This phenomenon is caused by the fact that the knowledge about the PHD (related to the posterior probability density) is gradually accumulated over time, and the uncertainty about the mean of mixture Gaussian distribution of PHD gradually decreases at the same time. This result proves that the proposed data assimilation algorithm can make the uncertainty of the estimated results be controlled in a limited range and it can overcome the uncertainty of the measurements including clutters, false alarms, missed detections, and noises.Figure 10
PHD of single pirate agent versus time: (a) surface plot of the distribution of PHD; (b) contours of the marginal PHD.
(a)
(b)To capture the average performance, we run 1000 Monte Carlo trials for each experiment with the same setup. Figure11 shows the Monte Carlo average of the mean squared error at every step for the assimilated pirate states and the simulated ones. We can find that the mean squared errors between the real pirate states and the assimilated ones are smaller than those of the simulated ones. As time evolves, the mean squared errors between the real pirate states and the simulated ones gradually increase, while the mean squared errors between the real ones and the assimilated ones are stable and limited. This figure reflects the differences between using data assimilation and not using it.Figure 11
1000 MC run average mean squared errors versus time for assimilated and simulated results.At each time step, the data assimilation algorithm assimilated the measurements by using the RFS based Gaussian mixture. To see the data assimilation process more clearly, we give the assimilation result at a single time step. Here we analyze the result at time step 57. To see the distribution of the PHD more intuitively, we give the surface plot of the distribution of PHD. From Figure12, we can know that the assimilated states of the pirate vessel are the states having the largest PHD. Since there is only one pirate vessel in the experiment, there is one two-dimensional Gaussian distribution in the xy plane.Figure 12
Surface plot of the distribution of assimilated PHD at time step 57.Figure13 shows the real pirate vessel state with cluttered measurements and the simulated one and the assimilated one at time step 57. From this figure, we can find that the assimilated state is much closer to the real one than the simulated one even though the measurements are heavily influenced by the clutters and noises. We also plot the mean value of the measurements. We can find that neither the mean value nor the measurements are closer to the real state than the assimilated one. This proves that the direct using of measurements without data assimilation is not applicable.Figure 13
The real state, assimilated state, simulated state, and corresponding measurements at time step 57.


### 6.3. Experiments on Multiple Pirate Vessels

Now consider a nonlinear bearing and range scenario with a total of 5 pirate vessels and a mothership. The birth process follows a Poisson RFS with PHD:(26)γkx=∑i=14wγiNx;mγi,Pγ,where wγ1 = wγ2 = 0.03 and wγ3 = wγ4 = 0.02, mγ1 = -1000,0,1500,0;0T, mγ2=-1000,0,-1000,0;0T, mγ3=1000,0,-1500,0;0T, mγ4=1000,0,500,0;0T, and Pγ=diag([50,50,50,50;6(π/180)]T)2.Like in Section6.2, we carry out three experiments with the same measurement model and simulation model. The configurations of the three experiments are identical except the turn rate. The real turn rates and erroneous ones of all 5 pirate agents and a mother ship are shown in Table 1. To show the ability of assimilating varying agents’ number of the proposed data assimilation algorithm, we make the pirate vessel agents’ number vary by setting different start and stop time.Table 1
Configurations of different pirate vessels.
Pirate vessels
Real turn rate (rad/s)
Biased turn rate (rad/s)
Start time (s)
Stop time (s)
Pirate 1
0.0058
0.0092
1
100
Pirate 2
0.0058
0.0092
1
80
Pirate 3
0.0087
0.0137
1
90
Pirate 4
0.0087
0.0137
30
100
Pirate 5
0.0058
0.0092
35
100
Mothership
0.0044
0.0069
1
50We assume that the merchant vessel is stationary at a particular location as a center. Thus, we can get the states in the Cartesian coordinates with the location of merchant vessels. Pirate vessels’ states are shown in Figure14.Figure 14
States of identical-twin experiments on multiple pirates.From Figure14, it is not difficult to find that, by assimilating the measurements, the assimilated pirate vessel states are closer to the real ones than the simulated ones.Figure15(a) shows the PHD of positions and Figure 15(b) shows the contours of the marginal PHD in all time steps. It can be learned that the PHD of all agents is gradually improved as the simulation time evolves. We can learn that, in the initial phase, the uncertainty of the PHD is large due to the limited and imprecise prior information, and as the simulation time evolves, the uncertainty of PHD gradually decreases. This phenomenon is consistent with the experiment for the single pirate agent. From this result, we can know that the varying agents’ number will not affect the ability of the proposed algorithm to control the uncertainty of estimated results in a limited range and overcome the uncertainty of measurements.Figure 15
PHD of multiple pirate agents versus time: (a) surface plot of the distribution of PHD; (b) contours of the marginal PHD.
(a)
(b)Unlike the analysis of single pirate vessel, the analysis of the multiple pirate vessels cannot use the distance error directly. Here we use the new metric which is proposed by Schuhmacher et al. for sets, named as the optimal subpattern assignment (OSPA) metric [45]. The OSPA metric is employed extensively as a performance evaluation criterion here.To capture the average performance, we also run 1000 Monte Carlo trials for each experiment with the same setup. Figure16 shows the Monte Carlo average of the OSPA distance, and we can find that the Monte Carlo average of the OSPA distance of simulated states gradually increases, while the OSPA distance of assimilated states is limited to a smaller range. This result reflects the fact that the multiple agent simulation using data assimilation is far more accurate than those not using it.Figure 16
1000 MC run average OSPA error versus time for simulated states and assimilated states.Figure17 shows the agents’ number in the real experiment and the assimilated agents’ number in the assimilated experiment. As the assimilated agents’ number is estimated by rounding off the mean of updated agents’ cardinality, a few estimation errors are caused. But mostly the assimilated numbers are consistent with the real agents’ numbers. This figure proves that the proposed data assimilation algorithm not only can assimilate the simulation states of the agents but also can estimate the varying agents’ number in the simulation. It is difficult for the standard PF based data assimilation algorithms to estimate the agents’ number at every time step.Figure 17
True versus assimilated cardinality over time.To describe the working mechanism of the proposed data assimilation algorithm more clearly, we give the assimilated results at a single time step. Here, we analyze time step 48 at which all the six pirate vessels are concurrent. To see the distribution of the PHD more intuitively, we give the surface plot of the distribution of PHD. From Figure18, we can know that the assimilated states of the pirate vessels are those with the largest PHD, and the number of the two-dimensional Gaussian distributions equals the number of pirate vessels.Figure 18
Surface plot of the distribution of resulting PHD at time step 48.Figure19 shows the real pirate vessels’ states, the simulated ones, the assimilated ones, the clutters, and noises. We can see that the assimilated pirate vessels’ states are closer to the real ones than the simulated ones. It can be seen that the assimilated states are not affected by the disorganized distribution of measurements.Figure 19
The real states, assimilated states, simulated states, and corresponding measurements at time step 48.


## 6.1. Experimental Design

We take advantage of the identical-twin experiment introduced in [31, 32] to verify the proposed data assimilation algorithm in ideal situations. In the identical-twin experiment, the simulation whose results are considered as real states is first executed with corresponding measurements which are collected during the process. And the measurements here are treated as the real measurements, since they are collected from the real states. Accordingly, we assimilate the measurements into DDDS by using the proposed data assimilation algorithm and then check whether the assimilated results are close to the real system states or not.We utilize three terms in the experiments, real pirate movement, assimilated pirate movement, and simulated pirate movement, to represent the results of the experiments. A real pirate movement is the simulated pirate movement from which the online measurements are obtained. A simulated pirate movement is the simulation states based on some erroneous parameters, for example, imprecise turn rate, which aims to comply with the fact that pirate movement simulation usually relies on imprecise parameters compared to real pirate movement. Finally, an assimilated pirate movement is the data assimilation enhanced simulation states based on the same erroneous parameters as in the simulated pirate movement. The goal of these experiments is to verify that an assimilated pirate movement creates more accurate simulation states by assimilating online measurements from the real pirate movement, even if it uses the erroneous parameters in the simulated pirate movement.


## 6.2. Experiments on Single Pirate Vessel in Clutter

The state space for the pirate vessel is four-dimensional with two position dimensions in the normalx,y Cartesian coordinate system and two velocity coordinates, x˙,y˙. We establish the 2D horizontal turn motion model relying on ship kinematics. We use a widely used ship motion model [44]. The dynamics equation is as follows:(23)xk=xk-1+Tvk-1sin⁡ϕk-1,yk=yk-1+Tvk-1cos⁡ϕk-1,ϕk=ϕk-1+Tvk-1Ωk-1+12Ωk-1-Ω0Tvk-1τeTvk-1τ,Ωk=Ωk-1eTvk-1τ+Ω01-eTvk-1τ,vk-1=K·v0=v01+1.9Ωk-12L2-1;here, τ=-p/2+p2/4-q/L, Ω0=Ω/v0, and x,y, Ω, and ϕ are ship position, velocity vector turn rate, and heading, respectively. T is the sampling interval.The simulation modelxk=φkxk-1,vk-1 of the vessel motion is nonlinear, where φk is the nonlinear state transition functions and vk-1 is an independent zero-mean Gaussian process noise; and Qk-1 is the covariance matrix. The vessel simulation model can be reformulated as(24)xk=φk-1xk-1+G·ϕ˙k-1,ϕ˙k=ϕ˙k-1+T·uk-1,where G=T2/2T0000T2/2TT,ϕ˙k-1~N·;0,σϕ˙2I,uk-1~N·;0;σu2I, T = 1 s, σϕ˙ = 15 m/s2, and σu=π/180 rad/s. The probability of survival is pS,kx = 0.99. The probability of detection is pD,kx = 0.98. The speed of the pirate vessel is 25–30 knots.The sensor measurement model in the experiments is as follows:(25)zk=arctan⁡xkykxk2+yk2+εk,where measurement noise εk~N·;0,Rk with Rk=diag(σθ2,σr2T); in the simulation experiments, we assume that σθ = 2(π/180) rad, and σr = 10 m.The false detection obeys the Poisson RFS with PHDλc = 3.2 × 10−3 (rad·m)−1 over the disc [0, 2π] rad × [0,2000] m.In our experiments, the differences between a real pirate movement and a simulated one are caused by the imprecise turn rateϕ˙. We use the imprecise turn rate as the erroneous parameters. The real turn rate is π/60≈0.0524 rad/s with random variances added every second. Our first experiment case uses the real turn rate to run the agent based simulation system and collects the measurements at every time step. Our second experiment case introduces errors to the turn rate and makes the other configurations exactly the same as the first experiment. The erroneous turn rate is 0.0524+0.01/3≈0.05573 rad/s with the same random variances in the first experiment. Our third experiment case uses the erroneous turn rate as the second experiment and the other configurations remain the same as the other two experiments. In this experiment case, the measurements collected in the first experiment are used as well as the proposed data assimilation algorithm.We assume that the merchant vessel is stationary at a particular location, which is the center in the Cartesian coordinates. In this way, we can get the pirate’s states. The states of the pirate vessel generated by the three experiments are shown in Figure9.Figure 9
States of identical-twin experiments on single pirate vessel in clutter.This set of experiments verifies the effectiveness of data assimilation algorithm when the turn rate has errors. From Figure9, we know that the real pirate vessel movement and the simulated one have large deviations due to the turn rate errors. By assimilating the measurements, the assimilated pirate vessel movement is obtained. From this figure, it is observed that the assimilated pirate vessel movement is much closer to the real one than the simulated one. Specifically, due to the erroneous turn rate, the simulated pirate states gradually deviate from the real pirate states. Using the proposed data assimilation algorithm, the assimilated pirate states overcome this problem and match the real pirate states with minor differences.Figure10(a) shows the PHD of positions and Figure 10(b) shows the contours of the marginal PHD in all time steps. We can learn that the PHD is gradually improved as the simulation time evolves. We can find that, in the initial phase, the uncertainty of the PHD is large due to the limited and imprecise prior information, and as the simulation time evolves, the uncertainty gradually decreases. This phenomenon is caused by the fact that the knowledge about the PHD (related to the posterior probability density) is gradually accumulated over time, and the uncertainty about the mean of mixture Gaussian distribution of PHD gradually decreases at the same time. This result proves that the proposed data assimilation algorithm can make the uncertainty of the estimated results be controlled in a limited range and it can overcome the uncertainty of the measurements including clutters, false alarms, missed detections, and noises.Figure 10
PHD of single pirate agent versus time: (a) surface plot of the distribution of PHD; (b) contours of the marginal PHD.
(a)
(b)To capture the average performance, we run 1000 Monte Carlo trials for each experiment with the same setup. Figure11 shows the Monte Carlo average of the mean squared error at every step for the assimilated pirate states and the simulated ones. We can find that the mean squared errors between the real pirate states and the assimilated ones are smaller than those of the simulated ones. As time evolves, the mean squared errors between the real pirate states and the simulated ones gradually increase, while the mean squared errors between the real ones and the assimilated ones are stable and limited. This figure reflects the differences between using data assimilation and not using it.Figure 11
1000 MC run average mean squared errors versus time for assimilated and simulated results.At each time step, the data assimilation algorithm assimilated the measurements by using the RFS based Gaussian mixture. To see the data assimilation process more clearly, we give the assimilation result at a single time step. Here we analyze the result at time step 57. To see the distribution of the PHD more intuitively, we give the surface plot of the distribution of PHD. From Figure12, we can know that the assimilated states of the pirate vessel are the states having the largest PHD. Since there is only one pirate vessel in the experiment, there is one two-dimensional Gaussian distribution in the xy plane.Figure 12
Surface plot of the distribution of assimilated PHD at time step 57.Figure13 shows the real pirate vessel state with cluttered measurements and the simulated one and the assimilated one at time step 57. From this figure, we can find that the assimilated state is much closer to the real one than the simulated one even though the measurements are heavily influenced by the clutters and noises. We also plot the mean value of the measurements. We can find that neither the mean value nor the measurements are closer to the real state than the assimilated one. This proves that the direct using of measurements without data assimilation is not applicable.Figure 13
The real state, assimilated state, simulated state, and corresponding measurements at time step 57.


## 6.3. Experiments on Multiple Pirate Vessels

Now consider a nonlinear bearing and range scenario with a total of 5 pirate vessels and a mothership. The birth process follows a Poisson RFS with PHD:(26)γkx=∑i=14wγiNx;mγi,Pγ,where wγ1 = wγ2 = 0.03 and wγ3 = wγ4 = 0.02, mγ1 = -1000,0,1500,0;0T, mγ2=-1000,0,-1000,0;0T, mγ3=1000,0,-1500,0;0T, mγ4=1000,0,500,0;0T, and Pγ=diag([50,50,50,50;6(π/180)]T)2.Like in Section6.2, we carry out three experiments with the same measurement model and simulation model. The configurations of the three experiments are identical except the turn rate. The real turn rates and erroneous ones of all 5 pirate agents and a mother ship are shown in Table 1. To show the ability of assimilating varying agents’ number of the proposed data assimilation algorithm, we make the pirate vessel agents’ number vary by setting different start and stop time.Table 1
Configurations of different pirate vessels.
Pirate vessels
Real turn rate (rad/s)
Biased turn rate (rad/s)
Start time (s)
Stop time (s)
Pirate 1
0.0058
0.0092
1
100
Pirate 2
0.0058
0.0092
1
80
Pirate 3
0.0087
0.0137
1
90
Pirate 4
0.0087
0.0137
30
100
Pirate 5
0.0058
0.0092
35
100
Mothership
0.0044
0.0069
1
50We assume that the merchant vessel is stationary at a particular location as a center. Thus, we can get the states in the Cartesian coordinates with the location of merchant vessels. Pirate vessels’ states are shown in Figure14.Figure 14
States of identical-twin experiments on multiple pirates.From Figure14, it is not difficult to find that, by assimilating the measurements, the assimilated pirate vessel states are closer to the real ones than the simulated ones.Figure15(a) shows the PHD of positions and Figure 15(b) shows the contours of the marginal PHD in all time steps. It can be learned that the PHD of all agents is gradually improved as the simulation time evolves. We can learn that, in the initial phase, the uncertainty of the PHD is large due to the limited and imprecise prior information, and as the simulation time evolves, the uncertainty of PHD gradually decreases. This phenomenon is consistent with the experiment for the single pirate agent. From this result, we can know that the varying agents’ number will not affect the ability of the proposed algorithm to control the uncertainty of estimated results in a limited range and overcome the uncertainty of measurements.Figure 15
PHD of multiple pirate agents versus time: (a) surface plot of the distribution of PHD; (b) contours of the marginal PHD.
(a)
(b)Unlike the analysis of single pirate vessel, the analysis of the multiple pirate vessels cannot use the distance error directly. Here we use the new metric which is proposed by Schuhmacher et al. for sets, named as the optimal subpattern assignment (OSPA) metric [45]. The OSPA metric is employed extensively as a performance evaluation criterion here.To capture the average performance, we also run 1000 Monte Carlo trials for each experiment with the same setup. Figure16 shows the Monte Carlo average of the OSPA distance, and we can find that the Monte Carlo average of the OSPA distance of simulated states gradually increases, while the OSPA distance of assimilated states is limited to a smaller range. This result reflects the fact that the multiple agent simulation using data assimilation is far more accurate than those not using it.Figure 16
1000 MC run average OSPA error versus time for simulated states and assimilated states.Figure17 shows the agents’ number in the real experiment and the assimilated agents’ number in the assimilated experiment. As the assimilated agents’ number is estimated by rounding off the mean of updated agents’ cardinality, a few estimation errors are caused. But mostly the assimilated numbers are consistent with the real agents’ numbers. This figure proves that the proposed data assimilation algorithm not only can assimilate the simulation states of the agents but also can estimate the varying agents’ number in the simulation. It is difficult for the standard PF based data assimilation algorithms to estimate the agents’ number at every time step.Figure 17
True versus assimilated cardinality over time.To describe the working mechanism of the proposed data assimilation algorithm more clearly, we give the assimilated results at a single time step. Here, we analyze time step 48 at which all the six pirate vessels are concurrent. To see the distribution of the PHD more intuitively, we give the surface plot of the distribution of PHD. From Figure18, we can know that the assimilated states of the pirate vessels are those with the largest PHD, and the number of the two-dimensional Gaussian distributions equals the number of pirate vessels.Figure 18
Surface plot of the distribution of resulting PHD at time step 48.Figure19 shows the real pirate vessels’ states, the simulated ones, the assimilated ones, the clutters, and noises. We can see that the assimilated pirate vessels’ states are closer to the real ones than the simulated ones. It can be seen that the assimilated states are not affected by the disorganized distribution of measurements.Figure 19
The real states, assimilated states, simulated states, and corresponding measurements at time step 48.


## 7. Conclusions

In this paper, a RFS based data assimilation algorithm is proposed for the agent based DDDS of maritime pirate activity. The main limitations of the traditional PF based and Kalman filter based data assimilation algorithms have been analyzed; therefore by employing RFS these limitations can be effectively overcome. This algorithm bears the ability to meet the demands of practical applications. Additionally, we designed the RFS based measurement model and simulation model for the maritime sensors and vessel agents. We implemented the RFS based data assimilation algorithm by using RFS based Gaussian mixture. We tested the effectiveness of the proposed algorithm by using the identical-twin experiments of maritime pirate activity. The experiment results prove that the proposed data assimilation algorithm can assimilate the input measurements and improve the accuracy of simulation results to a certain extent.



---

*Source: 5307219-2017-12-11.xml*