Nuclear Accident Emergency Response System: Radiation Field Estimation and Evacuation

Abstract

In this paper, a nuclear accident emergency response system based on unmanned aerial vehicles (UAVs) and bus collaboration is designed for radiation field estimation and evacuation. When a nuclear accident occurs, the radiation field is estimated firstly using the measurements acquired by UAVs. Based on the Cramer–Rao Lower Bound (CRLB), the coordinate optimization combined with UAV routing is formulated as a nonconvex mixed integer nonlinear programming (MINLP) problem to maximize the estimation accuracy. Further, a two-stage solution procedure based on genetic algorithm (GA) is proposed to solve the above problem. Then, taking the predicted radiation field as input, personnel in the emergency planning zone (EPZ) are evacuated to shelters by buses. The evacuation routing problem for minimizing both the total radiation exposure and evacuation time is formulated as a mixed integer linear programming (MILP) problem, which is solvable with efficient commercial solvers, such as Gurobi and CPLEX. The simulation results indicate that the system can provide effective help for emergency management under the nuclear accident scenarios.

1. Introduction

The nuclear accidents at Chernobyl in 1986 and Fukushima Daiichi in 2011 have caused serious casualties, economic losses, and environmental impacts [1,2,3]. Therefore, the emergency response to nuclear accidents, as the last barrier of nuclear safety, seems to be extremely important for the sustainable development of ecology and economy. With the consideration of limited resources, constrained time and radiation exposure risks to evacuees in the nuclear accident, the emergency response work is mainly divided into two parts, one is to estimate the radiation field for the risk assessment, and the other is to evacuate residents to shelters as soon as possible.

The radiation field estimation mainly uses the atmospheric diffusion simulation. However, the simulation accuracy depends heavily on the fidelity of model and the parameters. Due to the dynamic nature of the atmospheric diffusion, the parameters cannot be obtained accurately, resulting in lower estimation accuracy of the diffusion model. Based on real-time measurement, the estimation accuracy is improved by updating parameters in real time [4]. To obtain valid real-time measurement, UAVs that provide mobility for sensors have been widely used in air contaminant monitoring [5,6,7].

Based on sufficient and accurate information, including radiation field, population distribution, and so on, a reasonable evacuation plan can be formulated to ensure the safety of evacuees. In order to reduce the evacuation time and casualties, various methods, such as contraflow [8], traffic pattern simulation [9,10], and pedestrian dynamics [11], have been proposed, in which personnel are evacuated by personal vehicles. However, due to the large-scale population and limited road capacities in the area of the nuclear accident, bus-based evacuations are frequently used in emergency evacuation [12,13].

In this paper, a nuclear accident emergency response system based on UAV and bus collaboration is designed. As shown in Figure 1, the framework of this system is composed of radiation field estimation and emergency evacuation. Firstly, UAVs carry sensors to measure the radiation concentrations at the optimal coordinates that contain the maximum amount of information to maximize the estimation accuracy of radiation field. In detail, the amount of information is quantified by a CRLB-base metric and the coordinate optimization combined with UAV routing is formulated as a noncovex MINLP problem. A set of binary variables, specifically, the division of the measurement coordinate set for each UAV, ensure the coordinates accessibility with consideration of the speed of UAV. Then, the MINLP problem is solved by a two-stage solution procedure based on GA [14]. According to the measurements obtained by UAVs, maximum likelihood estimate (MLE) [15] is applied to update the parameters of diffusion model describing the radiation field. This process is repeated for prescribed rounds to accurately estimate the time-varying parameters. Taking the predicted radiation field as input, the bus-based evacuation planning is formulated as an MILP problem to minimize both the evacuation time and the radiation exposure to evacuees, and then the MILP problem is directly solved by the commercial solver such as Gurobi and CPLEX. Finally, the buses carry out the evacuation plan.

There are three key contributions in this paper:

1.

Design a nuclear accident emergency response system consisting of radiation field estimation and evacuation. Based on UAV and bus collaboration, it can provide efficient, reliable and safe nuclear emergency response strategy for the decision-maker;

2.

Analyze the optimal measurement coordinates considering the mobility of UAV. The CRLB-based coordinates optimization combined with UAV routing is formulated as an MINLP problem, and then solved by a two-stage solution procedure;

3.

Improve the bus evacuation MILP model proposed by Bolia [16], in which both the evacuation time and the radiation exposure to evacuees are taken into consideration. The optimal evacuation route for buses can be directly obtained by commercial solvers.

The remainder of this paper is organized as follows. Section 2 lays out the literature review. In Section 3, the radiation measurement model and the CRLB-based metric of parameter estimation accuracy are introduced. Then, the coordinates optimization problem with UAV routing is formulated and a two-stage solution procedure based on GA is proposed. In Section 4, some assumptions and descriptions are given firstly, and then the mathematical formulation for the bus evacuation problem is presented. Simulation results illustrate the effectiveness of the proposed nuclear accident emergency response system in Section 5. Section 6 summarizes this paper.

2. Literature Review

In terms of atmospheric diffusion simulation, various diffusion models have been proposed and used in the prediction of air contaminant dispersion, which describe the physical process of atmospheric diffusion using a mathematical expression. These models include Gaussian model [17,18], Lagrangian model [19], and computational fluid dynamics [20]. Wang et al. [21] propose a real-time data driven simulation of Gaussian puff model based on UAV sensory system. Hiemstra et al. [22] utilize particle filter to assimilate the observations of radiation into an atmospheric transport model for more accurate simulation. Fang et al. [23] compare the maximum likelihood-expectation maximization and algebraic reconstruction algorithms for reconstruction quality of Gaussian plume model. It can be noticed that the structure of diffusion model has been widely researched [24]. However, the parameters that are always regarded as a prior knowledge are also the decisive factor for simulation accuracy, and the parameters estimation has not been well studied.

As is well known, the CRLB [25] is widely used as the theoretical metric that quantifies the accuracy of parameter estimation. Ristic et al. [26] present theoretical analysis of the achievable accuracy in estimation of the radiation source term parameters including location and release amount. In addition, the framework of CRLB method can also be applied to the target localization problem [27,28,29,30,31]. For instance, Yang et al. [27] propose a Fisher information matrix (FIM)-based metric to optimize coordination strategy for target tracking under additive and multiplicative noises. Xu et al. [28] develop two new closed-form solutions to minimize the trace of CRLB in order to achieve optimal localization performance.

Once the radiation field is obtained, personnel in the EPZ need to be evacuated to shelters as soon as possible. There are extensive investigations on evacuation planning in various man-made or natural disasters, such as hurricanes, earthquakes, flooding, landslide, and nuclear accidents [32,33,34,35,36,37,38]. Dikas et al. [35] propose a new two-index formulation for the bus evacuation problem and its variants. What is more, a hybrid solution framework that integrated of large neighborhood search (LNS), variable neighborhood search (VNS), and Column Generation is designed. He et al. [37] adopt Benders decomposition to address the dynamic resource allocation problem for evacuation planning of large-scale networks. Goerigk et al. [38] take the individual traffic, public traffic and location decision into consideration, and propose a genetic algorithm to solve the evacuation problem.

In order to ensure the safety of the evacuees under the nuclear accident scenarios, the radiation exposure should be considered into the evacuation problem. Urbanik et al. [39] conduct the evacuation time estimate analysis to provide information relevant to the development of effective evacuation plans. Tan et al. [40] design a nuclear emergency parallel evacuation system, in which the maximum evacuation time for each demand is determined by the radiation exposure risk. Huang et al. [41] develop an inexact fuzzy stochastic chance constrained programming method to address the radiation uncertainty in the nuclear accident. However, the evacuation problem considering the maximum radiation dose for evacuees has not been well studied [42].

3. UAV-Based Nuclear Radiation Field Estimation

This section solves the problem of UAV-based nuclear radiation field estimation. We firstly introduce the radiation measurement model with white Gaussian noise. Then, a CRLB-based metric, which depends on the measurement coordinate, is introduced to measure estimation accuracy. The coordinates optimization combined with UAV routing is formulated as an MINLP problem. Finally, a two-stage solution procedure based on GA is proposed to obtain the optimal feasible solution.

3.1. Radiation Measurement Model

When a nuclear accident occurs, the distribution of radiation must be time-varying. Being different from the stable diffusion model [23], a Gaussian puff model, which describes the instantaneous concentration of radiation, is adopted in this paper to preserve the time characteristic. The radiation source is located at the coordinate system’s origin. Then, the radiation concentration at a space-time coordinate $(x,y,z,t)$ is given by

$$C=\frac{2Q_0}{{\left(2\pi \right)}^{3/2}{\delta }_x{\delta }_y{\delta }_z}exp\left[-\frac{{\left(x-ut\right)}^2}{2{\delta }_x^2}-\frac{y^2}{2{\delta }_y^2}-\frac{z^2}{2{\delta }_z^2}\right],$$

(1)

where $exp(·)$ denotes the exponential function; $Q_0$ known as a prior denotes the total amount of leaked radioactive material; u denotes the constant wind speed, and the wind direction is along the X-axis; ${\delta }_x$, ${\delta }_y$, and ${\delta }_z$ are the diffusion parameters along the X-, Y-, and Z-axes.

After determining the structure of diffusion model, the parameters of them affect the estimation accuracy of radiation field. Although the diffusion parameters are always from the empirical formula in other papers [40], we consider them as parameters that need to be estimated accurately. In order to ensure the uniqueness of solution, the total number of measurements N is at least equal to the number of parameters to be estimated.

Assume that a group of UAVs carry sensors to measure the radiation concentrations in the radiation field. The coordinates of the i-th measurement $c_i$ are $(x_i,y_i,z_i,t_i)$. Due to the huge order of magnitude difference of radiation concentrations measured at different coordinates, it is assumed that there is a multiplicative error ${\nu }_i$ in measurement, i.e., $c_i=C_i(1+{\nu }_i)$. The i-th logarithmic measurement is modelled by

$$m_i=ln{c}_i=ln{C}_i+{\omega }_i,$$

(2)

where ${\omega }_i=ln(1+{\nu }_i)$ is additive white Gaussian noise with zero mean and covariance ${\sigma }^2$, i.e., ${\omega }_i\sim \mathcal{N}(0,{\sigma }^2)$. As a result, the i-th logarithmic measurement also satisfies Gaussian distribution $m_i\sim \mathcal{N}(ln{C}_i,{\sigma }^2)$.

The parameter vector $\theta$ is denoted as $({\xi }_1,{\xi }_2,{\xi }_3)$, which equals to $(\frac{1}{{\delta }_x},\frac{1}{{\delta }_y},\frac{1}{{\delta }_z})$. According to (2), the probability density function (PDF) of $m_i$ is obtained as

$$p\left(m_i\mid \mathbf{\theta }\right)=\frac{1}{\sqrt{2\pi }\sigma }exp\left[-\frac{{\left(m_i-ln{C}_i\left(\theta \right)\right)}^2}{2{\sigma }^2}\right].$$

(3)

Due to the independence of measurements, the vector consisting of N measurements follows N-dimensional Gaussian distribution, i.e.,

$$\mathcal{M}\sim \mathcal{N}(\mathsf{\Phi },\mathsf{\Sigma }),$$

(4)

where $\mathcal{M}={\left[m_1,m_2,\dots ,m_N\right]}^{\top }$, $\mathsf{\Phi }=[ln{C}_1,ln{C}_2,\dots ,$$ln{C}_N{]}^{\top }$, and $\mathsf{\Sigma }=diag{\left\{{\sigma }^2,{\sigma }^2,\dots ,{\sigma }^2\right\}}_{N\times 1}$. Thus, the joint PDF can be derived from (3) as

$$p\left(\mathcal{M}\mid \mathbf{\theta }\right)=\prod _{i=1}^{N}p\left(m_i\mid \mathbf{\theta }\right)=\frac{1}{{\left(2\pi {\sigma }^2\right)}^{\frac{N}{2}}}exp\left[-\frac{1}{2{\sigma }^2}\sum _{i=1}^N{\left(m_i-ln{C}_i\left(\theta \right)\right)}^2\right].$$

(5)

3.2. CRLB-Based Metric

As is well known, the FIM stands for the information available in measurements [27]. A larger FIM means more information can be retrieved from measurements, thus, the accuracy of parameter estimation is improved. Let $\mathit{J}$ denote the FIM, and from the definition in [25], we have

$$\mathit{J}=-\mathrm{E}\left\{{\nabla }_{\theta }{\nabla }_{\theta }^{\top }lnp(M\mid \theta )\right\},$$

(6)

where $\mathrm{E}[·]$ is the mathematical expectation operator; ${\nabla }_{\theta }$ is gradient operator with respect to $\theta$; Operator ${\nabla }_{\theta }{\nabla }_{\theta }^{\top }\equiv {\Delta }_{\theta }$ is the Hessian matrix with respect to $\theta$.

For the logarithmic measurement $m_i,i=1,\dots ,N$, the general expression of $\mathit{J}$ [25] can be written as

$$\mathit{J}=\frac{1}{{\sigma }^2}\sum _{i=1}^N\left[{\nabla }_{\theta }ln{C}_i\left(\theta \right)\right]·[{\nabla }_{\theta }^{\top }ln{C}_i\left(\theta \right)],$$

(7)

where the complete expressions for the components of the vector ${\nabla }_{\theta }ln{C}_i\left(\theta \right)$ are given by

$$\begin{array}{c}\hfill \partial ln{C}_i\left(\theta \right)/\partial {\xi }_1=\frac{1}{{\xi }_1}-{(x_i-u{t}_i)}^2{\xi }_1,\\ \hfill \partial ln{C}_i\left(\theta \right)/\partial {\xi }_2=\frac{1}{{\xi }_2}-y_i^2{\xi }_2,\\ \hfill \partial ln{C}_i\left(\theta \right)/\partial {\xi }_3=\frac{1}{{\xi }_3}-z_i^2{\xi }_3.\end{array}$$

(8)

We define

$$\begin{array}{c}\hfill \mathit{\alpha }=\frac{1}{\sigma }{\left[\frac{1-{(x_1-u{t}_1)}^2{\xi }_1^2}{{\xi }_1},\dots ,\frac{1-{(x_N-u{t}_N)}^2{\xi }_1^2}{{\xi }_1}\right]}^{\top },\\ \hfill \mathit{\beta }=\frac{1}{\sigma }{\left[\frac{1}{{\xi }_2}-y_1^2{\xi }_2,\dots ,\frac{1}{{\xi }_2}-y_N^2{\xi }_2\right]}^{\top },\\ \hfill \mathit{\gamma }=\frac{1}{\sigma }{\left[\frac{1}{{\xi }_3}-z_1^2{\xi }_3,\dots ,\frac{1}{{\xi }_3}-z_N^2{\xi }_3\right]}^{\top }.\end{array}$$

(9)

Then, the FIM $\mathit{J}$ can be rewritten as

$$\mathit{J}=\left[\begin{array}{ccc}{\mathit{\alpha }}^{\top }\mathit{\alpha }& {\mathit{\alpha }}^{\top }\mathit{\beta }& {\mathit{\alpha }}^{\top }\mathit{\gamma }\\ {\mathit{\beta }}^{\top }\mathit{\alpha }& {\mathit{\beta }}^{\top }\mathit{\beta }& {\mathit{\beta }}^{\top }\mathit{\gamma }\\ {\mathit{\gamma }}^{\top }\mathit{\alpha }& {\mathit{\gamma }}^{\top }\mathit{\beta }& {\mathit{\gamma }}^{\top }\mathit{\gamma }\end{array}\right]$$

(10)

It is noted that $\mathit{J}$ as a matrix is not convenient for further analysis. Therefore, the trace of CRLB is used as the accuracy metric in this paper, which is in fact the well-known A-optimality criterion [43]. By minimizing the metric, the information retrieved from measurements is maximized.

Let ${\vartheta }_1$, ${\vartheta }_2$ and ${\vartheta }_3$ denote the intersection angle between vectors $\mathit{\alpha }$ and $\mathit{\beta }$, $\mathit{\alpha }$ and $\mathit{\gamma }$, and $\mathit{\beta }$ and $\mathit{\gamma }$. In addition, $\parallel ·\parallel$ denotes the norm of vector, and for example, ${\mathit{\alpha }}^{\top }\mathit{\beta }=\parallel \mathit{\alpha }\parallel \parallel \mathit{\beta }\parallel cos{\vartheta }_1$. According to (10) and [25], the CRLB can be calculated as

$$\mathit{CRLB}={\mathit{J}}^{-1},$$

(11)

As rigorously proved in [28], we have

$$tr\left(\mathit{CRLB}\right)\ge \frac{1}{{\parallel \mathit{\alpha }\parallel }^2}+\frac{1}{{\parallel \mathit{\beta }\parallel }^2}+\frac{1}{{\parallel \mathit{\gamma }\parallel }^2},$$

(12)

the above equality holds when

$$cos{\vartheta }_1=cos{\vartheta }_2=cos{\vartheta }_3=0,$$

(13)

which means the CRLB and FIM are both diagonal. Therefore, minimizing the trace of CRLB is equivalent to minimizing the the right term of (12) and subject to (13). Replacing $\mathit{\alpha }$, $\mathit{\beta }$, $\mathit{\gamma }$, ${\vartheta }_1$, ${\vartheta }_2$ and ${\vartheta }_3$ by $x_i,y_i,z_i,t_i$, the metric $CM$ can be explicitly given by

$$\begin{array}{cc}\hfill & CM=\frac{1}{{\parallel \mathit{\alpha }\parallel }^2}+\frac{1}{{\parallel \mathit{\beta }\parallel }^2}+\frac{1}{{\parallel \mathit{\gamma }\parallel }^2}\hfill \\ \hfill & =\frac{1}{{\sigma }^2}\left[{\left(\sum _{i=1}^N{(\frac{1}{{\xi }_1}-{(x_i-u{t}_i)}^2{\xi }_1)}^2\right)}^{-1}+{\left(\sum _{i=1}^N{(\frac{1}{{\xi }_2}-y_i^2{\xi }_2)}^2\right)}^{-1}+{\left(\sum _{i=1}^N{(\frac{1}{{\xi }_3}-z_i^2{\xi }_3)}^2\right)}^{-1}\right],\hfill \end{array}$$

(14)

with equality constraints,

$$\begin{array}{c}\hfill \sum _{i=1}^N\left(1-{(x_i-u{t}_i)}^2{\xi }_1^2\right)\left(1-y_i^2{\xi }_2^2\right)=0,\\ \hfill \sum _{i=1}^N\left(1-{(x_i-u{t}_i)}^2{\xi }_1^2\right)\left(1-y_i^2{\xi }_2^2\right)=0,\\ \hfill \sum _{i=1}^N\left(1-y_i^2{\xi }_2^2\right)\left(1-z_i^2{\xi }_3^2\right)=0.\end{array}$$

(15)

3.3. Coordinates Optimization Problem Formulation

According to the CRLB-based metric, the accuracy of parameter estimation depends on the measurement coordinates. Our objective is to optimize the measurement coordinates to maximize the estimation accuracy under the constraints on measurement number and mobility of UAV.

Taking into account the speed of UAV v, a set partition problem is integrated to ensure the accessibility of the optimal measurement coordinates. Let $\mathcal{S}$ be the set of UAV depots, $\mathcal{K}$ be the set of UAVs and $\mathcal{N}$ be the set of measurement coordinates. The binary variable $g_{ik}$, which is equal to 1 if a coordinate i is reached by UAV k and 0 otherwise, is introduced to divide $\mathcal{N}$ into small sets for each UAV.

In order to find a feasible solution, the equality constraints in (15) are relaxed into inequality constraints by setting a small threshold $ϵ$. Let ${\pi }_i$ denote the i-th measurement coordinate composed of decision variables $x_i,y_i,z_i,t_i$. The UAV-based coordinates optimization is formulated as an MINLP problem:

$$\underset{\{x_i,y_i,z_i,t_i\},i\in \mathcal{N}}{argmin}CM$$

(16)

s.t.

$$\left|\sum _{i\in \mathcal{N}}\left(1-{(x_i-u{t}_i)}^2{\xi }_1^2\right)\left(1-y_i^2{\xi }_2^2\right)\right|\le ϵ$$

(17)

$$\left|\sum _{i\in \mathcal{N}}\left(1-{(x_i-u{t}_i)}^2{\xi }_1^2\right)\left(1-z_i^2{\xi }_3^2\right)\right|\le ϵ$$

(18)

$$\left|\sum _{i\in \mathcal{N}}\left(1-y_i^2{\xi }_2^2\right)\left(1-z_i^2{\xi }_3^2\right)\right|\le ϵ$$

(19)

$$\sum _{k\in \mathcal{K}}{g}_{ik}=1,\forall i\in \mathcal{N}\cup \mathcal{S}$$

(20)

$$-M(1-g_{ik})-M(1-g_{jk})+\parallel {\pi }_i-{\pi }_j\parallel \le |t_j{g}_{jk}-t_i{g}_{ik}|v,\forall i,j\in \mathcal{N}\cup \mathcal{S},k\in \mathcal{K}$$

(21)

$$g_{ik}\in \{0,1\},\forall i\in \mathcal{N},k\in \mathcal{K}$$

(22)

$${\pi }_l\le {\pi }_i\le {\pi }_u,\forall i\in \mathcal{N}$$

(23)

Constraints (17)–(19) are the relaxed inequality constraints from (15), where $|·|$ denotes the absolute value of scalar. Constraints (20) denote that every coordinate can only be measured by one UAV. Constraints (21) ensure that if coordinate i and coordinate j are reached by the same UAV k, then the distance between the two coordinates must be less than the maximum flight distance of UAV within the time interval. Otherwise, the first two penalty terms also guarantee the inequality, where M is a very large positive number. Constraints (22) are the binary restriction on $g_{ik}$. Constraints (23) indicate that the measurement coordinates are bounded, where ${\pi }_l$ and ${\pi }_u$ are the lower bound and upper bound, respectively.

3.4. Two-Stage Solution Procedure

The nonconvex MINLP problem in Section 3.2 is complicated, which cannot be solved by commercial solvers, such as CPLEX and Gurobi. Therefore, a two-stage solution procedure is proposed. In the first stage, the genetic algorithm [14] is utilized to determine the binary variable $g_{ik}$ and obtain a initial solution for continuous variable ${\pi }_i$. In the second stage, the problem is transformed to a nonlinear programming (NLP) problem after fixing the binary variable, then the coordinates are further optimized based on the initial solution.

The genetic algorithm [44], which originates from the computer simulation research on the biological system, is a random search and optimization method that mimics the development of nature biological evolutionary mechanisms. Based on the population selection, crossover, and mutation, it can automatically acquire and accumulate knowledge about search space during searching, and then the optimal integer variable is obtained. However, the metaheuristic cannot ensure the whole optimization converge to global minimum. Based on the initial solution of continuous variable and the fixed integer variable obtained from GA, Global Search in MATLAB is adopted to improve the solution quality of resultant NLP problem. Specifically, taking the local minimum from initial point as a benchmark, the optimal local minimum obtained from random generated starting points is obtained. With the two-stage solution procedure described below by Algorithm 1, the optimal feasible solution of UAV-based coordinates optimization problem can be obtained within a very short time.

Algorithm 1 Two-stage solution procedure.

Require:$Q_0$; ${\pi }_l$; ${\pi }_u$; v; the initial model parameters $\widehat{\mathbf{\theta }}$; the initial UAV depots $s_0$  Ensure:${\pi }_i$; $g_{ik}$    1: Solve the MINLP by GA to obtain the optimal integer variable $g_{ik}$ and an initial solution of measurement coordinates ${\pi }_i$;   2: Fix the integer variable $g_{ik}$ to obtain the optimal feasible measurement coordinates ${\pi }_i$ by Global Search;

As explained in [15], a maximum likelihood estimator can be applied to estimate the diffusion parameters maximizing the joint PDF with N UAV measurements in (5). The MLE of the diffusion parameters is

$$({\overline{\delta }}_x,{\overline{\delta }}_y,{\overline{\delta }}_z)=\underset{\{{\delta }_x,{\delta }_y,{\delta }_z\}}{argmin}f({\delta }_x,{\delta }_y,{\delta }_z),$$

(24)

where $f({\delta }_x,{\delta }_y,{\delta }_z)=\frac{1}{2{\sigma }^2}{\sum }_{i=1}^N{(m_i-ln{C}_i)}^2$. Then, the genetic algorithm is also utilized to solve the aforementioned unconstrained nonlinear optimization problem in (24).

The detailed procedure of the nuclear radiation field estimation method is described in Algorithm 2. By discretizing the time-varying parameters of diffusion model, our proposed method consisting of the UAV-based coordinates optimization and the MLE-based parameter optimization can accurately estimate the dynamic nuclear radiation field, which is proven by the simulation results in Section 4.

Algorithm 2 UAV-based nuclear radiation field estimation.

Require:$Q_0$; ${\pi }_l$; ${\pi }_u$; v; $\widehat{\mathbf{\theta }}$; $s_0$; the time interval $\Delta t$  Ensure: The time-varying parameters of diffusion model   1: repeat   2:       Solve the coordinates optimization problem by Algorithm 1;   3:       Measure the radiation concentrations at the optimal coordinates by UAVs   4:       Estimate the diffusion parameters $\overline{\mathbf{\theta }}$ according to (24);   5:       Set the last location of UAV as new depot $s_0$ according to the flow of UAV;   6:       Set $\widehat{\mathbf{\theta }}=\overline{\mathbf{\theta }}$, $t=t+\Delta t$   7: until The maximum number of rounds for parameter estimation is reached.

4. Bus-Based Nuclear Emergency Evacuation

This section solves the problem of bus-based nuclear emergency evacuation planning. Based on the radiation field acquired in the above section and other basic information, the inputs required for developing the evacuation problem are satisfied. In order to determine the optimal bus operating strategies while guaranteeing the evacuees’ safety during the evacuation process, an MILP model inspired by the work of Bolia [16] is formulated in this paper.

4.1. Assumption and Description

Before explaining the details of the model, the following assumptions are made.

1.

People arrive at the nearest pickup point in advance to wait for evacuation, and the transfer time and radiation exposure are ignored;

2.

The loading and unloading time of buses for pickup points and shelters are ignored;

3.

The capacities of buses at different depots and the demands of different pickup points are known;

4.

The locations of the depots, pickup points, and shelters are known and the travel time between them are constant;

5.

Each pickup point has a particular shelter, i.e., a bus only takes evacuees from a pickup point to the assigned shelter during one trip, even if not fully loaded;

6.

The shelter can accommodate all evacuees from the corresponding pickup points;

7.

The radiation dose per second of each route and each pickup point are constant and known during one trip.

The evacuation network is modelled as a directed graph $\mathcal{G}(\mathcal{V},\mathcal{E})$, where $\mathcal{V}$ and $\mathcal{E}$ denote the set of nodes and arcs. $\mathcal{V}$ includes a set of depots $\mathcal{D}$ where buses are initially departed from, a set of pickup points $\mathcal{P}$ where people are waiting for evacuation and a set of shelters $\mathcal{F}$, i.e., $\mathcal{V}=\{\mathcal{D},\mathcal{P},\mathcal{F}\}$. Particularly, $f_p$ denotes the assigned shelter of pickup point $p\in \mathcal{P}$. $\mathcal{E}=\left\{(d,p)\right|d\in \mathcal{D},p\in \mathcal{P}\}\cup \left\{(p,f_p)\right|p\in \mathcal{P}\}\cup \left\{(f,p)\right|f\in \mathcal{F},p\in \mathcal{P}\}$ represents the arcs connecting nodes in $\mathcal{V}$. The round trips of buses with index are denoted as $\mathcal{T}=\{1,\dots ,t_{max}\}$, in which $t_{max}$ is limited by the available time. An example of the evacuation network is illustrated in Figure 2. The yellow nodes, red nodes, and green nodes represent depots, pickup points, and shelters. Additionally, the dashed lines are unidirectional, while the solid lines are bidirectional.

Let ${\mathcal{B}}_d$ be the set of buses at depot $d\in \mathcal{D}$, $B_d$ be the number of buses located at depot d and $q_{d,b}$ be the capacity of bus b of depot d. The number of evacuees waiting for evacuation at a pickup point $p\in \mathcal{P}$ is denoted by $D_p$. Let $t_{d,p}$ denote the travel time from depot d to pickup point p and $t_{p,f_p}$ the travel time from pickup point p to its assigned shelter $f_p$.

The binary decision variable $s_{d,b}^t$ takes 1 if t-th trip of bus b that depart from depot d happens, and 0 otherwise. Another binary variable $l_{d,b,p}^t$ takes 1 if a bus b runs to pickup point p for its t-th trip, and 0 otherwise. The loading and unloading time are denoted by $lt$ and $ut$.

The round trip time of a bus b of depot d for its t-th trip to pickup point p is composed of pickup time $P{T}_{d,b,p}^t$ and send time $S{T}_{d,b,p}^t$. For all trips T, the send time represents the travel time from pickup point p to t-th trip shelter, which is defined as

$$S{T}_{d,b,p}^t=t_{p,f_p}{l}_{d,b,p}^t,1\le t\le t_{max}.$$

(25)

However, the pickup time represents the travel time from depot d to pickup point p if $t=1$, and the travel time from the shelter reached on the $(t-1)$-th trip to pickup point p otherwise, which is defined as

$$P{T}_{d,b,p}^t=\left\{\begin{array}{cc}\hfill & t_{d,p}{l}_{d,b,p}^t,t=1\hfill \\ \hfill & \sum _{g\in P}{t}_{p,f_g}{l}_{d,b,g}^{t-1}{l}_{d,b,p}^t,2\le t\le t_{max}.\hfill \end{array}\right.$$

(26)

In order to deal with the binary bi-linear terms in (26), which may cause intractable computational complexity, the linearization method described in [16] is adopted to eliminate the nonlinearity. We introduce another binary variable $h_{d,b,p,g}^t=l_{d,b,g}^{t-1}{l}_{d,b,p}^t$ which takes 1 if a bus has run to pickup point g during the $(t-1)$th trip and it runs to pickup point p during t-th trip, and 0 otherwise.

As stated earlier, it is assumed that ${\eta }_p^t$ and ${\eta }_{p,f_p}^t$, which denote the radiation dose per second of pickup point p during the t-th trip and the radiation dose per second between p and $s_p$ during the t-th trip, can be obtained from the predicted radiation field and regarded as constant. The total radiation exposure to evacuees transported by bus b of depot d for the t-th trip is composed of waiting radiation $W{R}_{d,b,p}^t$ and evacuation radiation $E{R}_{d,b,p}^t$, which are given by

$$W{R}_{d,b,p}^t=\sum _{tt=1}^t\sum _k\left(\left(P{T}_{d,b,k}^{tt}+S{T}_{d,b,k}^{tt}+\left(lt+ut\right)l_{d,b,p}^t\right){\eta }_p^{tt}{q}_{d,b}\right)-S{T}_{d,b,p}^t{\eta }_p^t{q}_{d,b},$$

(27)

$$E{R}_{d,b,p}^t=S{T}_{d,b,p}^t{\eta }_{p,f_p}^t{q}_{d,b}.$$

(28)

4.2. Mathematical Formulation

In this subsection, the objective function and constraints of the bus-based evacuation model are explained in detail. Let $E{T}_{total}$ be the total evacuation time and $R_{total}$ be the total radiation exposure to all evacuees. The bus-based evacuation planning is formulated as an MILP problem:

$$\underset{\{s_{d,b}^t,l_{d,b,p}^t,h_{d,b,p,g}^t\}}{argmin}E{T}_{total}+\frac{1}{L}{R}_{total}$$

(29)

s.t.

$$\sum _b{s}_{d,b}^1\ge B_d,\forall d\in \mathcal{D}$$

(30)

$$s_{d,b}^t\ge s_{d,b}^{t+1},\forall d\in \mathcal{D},b\in {\mathcal{B}}_d,t\in \mathcal{T}$$

(31)

$$\sum _p{l}_{d,b,p}^t\le s_{d,b}^t,\forall d\in \mathcal{D},b\in {\mathcal{B}}_d,t\in \mathcal{T}$$

(32)

$$\sum _d\sum _b\sum _t{s}_{d,b}^t{q}_{d,b}\ge \sum _p{D}_p$$

(33)

$$\sum _d\sum _b\sum _t{l}_{d,b,p}^t{q}_{d,b}\ge D_p,\forall p\in \mathcal{P}$$

(34)

$$2h_{d,b,p,g}^t\le l_{d,b,p}^t+l_{d,b,g}^{t-1},\forall d\in \mathcal{D},b\in {\mathcal{B}}_d,t\in \mathcal{T},p,g\in \mathcal{P}$$

(35)

$$h_{d,b,p,g}^t\ge l_{d,b,p}^t+l_{d,b,g}^{t-1}-1,\forall d\in \mathcal{D},b\in {\mathcal{B}}_d,t\in \mathcal{T},p,g\in \mathcal{P}$$

(36)

$$R_{total}\ge \sum _d\sum _b\sum _p\sum _t\left(W{R}_{d,b,p}^t+E{R}_{d,b,p}^t\right)$$

(37)

$$E{T}_{total}\ge \sum _{t=1}^{t_{max}}\sum _p\left(P{T}_{d,b,p}^t+S{T}_{d,b,p}^t+(lt+ut)l_{d,b,p}^t\right),\forall d\in \mathcal{D},b\in {\mathcal{B}}_d$$

(38)

$$s_{d,b}^t,l_{d,b,p}^t\in \{0,1\},\forall d\in \mathcal{D},b\in {\mathcal{B}}_d,t\in \mathcal{T}$$

(39)

$$h_{d,b,p,g}^t\in \{0,1\},\forall d\in \mathcal{D},b\in {\mathcal{B}}_d,t\in \mathcal{T},p,g\in \mathcal{P}$$

(40)

The objective is to minimize both $E{T}_{total}$ and $R_{total}$ to determine the optimal bus operating strategies while guaranteeing the evacuees’ safety. Due to the order of magnitude difference in time and radiation, a large enough value L is introduced to ensure the balance between the first and second terms of (29). Constraints (30) ensure that every bus of depot d is used for evacuation. Constraints (31) ensure that a bus only start the subsequent trip after completing its current trip. Constraints (32) dictate that a bus only arrive at one pickup point during one trip if the trip in fact happens. Constraints (33) ensure that the number of people delivered by all buses during all trips is greater than the total number of evacuees waiting at the pickup points. Constraints (34) ensure that all evacuees waiting at each pickup point are picked up. Constraints (35) and (36) are for consistency in the definition of variable $h_{d,b,p,g}^t$. Constraint (37) denotes the total radiation exposure to evacuees across all trips of all buses. Constraints (38) denote the maximum evacuation time across all trips among all buses, i.e., the total evacuation time. Constraints (39) and (40) specify the decision variables $s_{d,b}^t$, $l_{d,b,p}^t$ and the auxiliary variable $h_{d,b,p,g}^t$ are binary. Combining all the constraints and the objective function, the MILP problem can be directly solved by commercial solvers, such as CPLEX and Gurubi.

5. Simulation Results

This section demonstrates the performance of the proposed nuclear emergency response system by the simulation results. It is assumed that our system is activated 600 s after the nuclear leakage. The speed of UAV is regarded as a constant 25 m/s. Due to the limitations of the number and speed of UAV, only three UAVs carry sensors to measure over a 1000 × 1000 × 500 m critical area to predict the whole radiation field. The total amount of leaked radioactive material is $Q_0=6.5\times {10}^{12}$, and the logarithmic measurement noise yields ${\omega }_i\sim \mathcal{N}(0,0.01)$. Note that the true value of parameter ${\mathbf{\theta }}^{\ast }$ is assumed to be clearly known to verify the performance of radiation field estimation,

$$\left\{\begin{array}{cc}\hfill {\delta }_x^{\ast }& =0.001{(t-600)}^2+40,\hfill \\ \hfill {\delta }_y^{\ast }& =0.0005{(t-1000)}^2+80,\hfill \\ \hfill {\delta }_z^{\ast }& =0.1t.\hfill \end{array}\right.$$

(41)

and the time interval is $\Delta t=20$ s. At each $[t,t+\Delta t]$, the proposed MINLP is solved by the two-stage solution procedure to find the optimal feasible solution.

Table 1. Optimal measurement coordinates for the first round of UAVs.

UAVs	x (m)	y (m)	z (m)	t (s)
1	173.027	244.802	26.545	608.697
	300.989	184.739	121.332	616.847
	303.188	180.700	125.570	618.166
	318.715	174.929	120.500	619.946
2	−109.613	−141.109	49.005	613.547
	−65.560	−156.113	63.143	616.615
	−39.174	−159.356	69.787	618.815
	−41.453	−158.657	70.184	619.958
3	299.851	247.711	4.064	615.587
	256.180	284.296	56.412	619.909

shows the optimal measurement coordinates for the first round of UAVs. It can be inferred from the table that the UAVs are reasonably routed to carry out their measurement tasks within the time window. Then, the optimal parameters are determined based on UAV measurements.

Three kinds of measurement coordinates strategies are conducted for parameter estimation: (1) The proposed UAV-based measurement strategy, which is composed of 10 optimal coordinates with the maximum amount of information; (2) 10 fixed measurement coordinates strategy; (3) 20 fixed measurement coordinates strategy. The parameter estimation performance is shown in Figure 3. In the left column, we compare the parameter estimation errors of the three strategies, in which the red lines are for the parameter estimation errors based on the proposed strategy, and, likewise, the blue and purple lines are for the errors based on 10 and 20 fixed coordinates strategy. In the right column, we compare the parameter change trends of three strategies with the true value within $[600,1200]$ s, in which the red lines represent for the true value of parameters, the blue, green, and purple lines are the estimated trends based on the proposed strategy, 10 and 20 fixed coordinates strategy. Obviously, the parameter estimation errors of the proposed strategy are smaller than the compared strategy. Therefore, the proposed strategy has better parameter estimation performance.

Since the parameter estimation problem solved by GA may converge to a local minimum, the parameter root mean square (RMS) errors for the first round are obtained by averaging over $s=10$ Monte Carlo runs, which is defined as

$$e=\sqrt{\frac{1}{s}\sum _{i=1}^s\left[{\left({\overline{\delta }}_{x,i}-{\delta }_x^{\ast }\right)}^2+{\left({\overline{\delta }}_{y,i}-{\delta }_y^{\ast }\right)}^2+{\left({\overline{\delta }}_{z,i}-{\delta }_z^{\ast }\right)}^2\right]}.$$

(42)

Figure 4 shows the corresponding parameter RMS errors for three strategies, displayed as a function of the number of Monte Carlo runs $s=1,\dots ,10$ used for averaging. The following observations can be made:

1.

More information can be obtained from more measurements, resulting in smaller parameter RMS error;

2.

Based on the optimal measurement coordinates, even smaller parameter RMS error can be achieved with fewer measurements.

In order to validate the performance for radiation field estimation, the ground is selected as a reference for predicting the radiation concentration at a certain moment. Figure 5 illustrates the radiation concentration error based on the optimal measurement coordinates and 20 fixed measurement coordinates. Comparing Figure 5a,b, it can be noticed that our proposed strategy predicts the radiation concentration more efficiently and accurately. The reason for the advantage is that the radiation concentration is extremely sensitive to the parameter errors. Therefore, the estimation accuracy is significantly improved while the parameter errors are only slightly decreased.

After a period of parameter estimation based on UAV measurements, the subsequent diffusion parameters can be fitted. Once the basic information, including the radiation field, population distribution, and node locations, is obtained, the bus-based nuclear emergency evacuation is initiated. Assume that the EPZ covers the area within a radius of 2500 m of the radiation source. There are 1500 people waiting for evacuation at different pickup points. Each pickup point is assigned to its nearest shelter. Due to the resource and time limitation, the maximum number of trips for each bus is set to be four. The speed of bus is a constant 20 m/s, the loading and unloading time are both constant 10 s. As stated earlier, the radiation dose per second for each route and pickup points is assumed to be constant during one trip. To do so, an invariable round trip time, which is the average of all possible round trip time, is prescribed. Taking the above round trip time and location of nodes as input, we can obtain ${\eta }_p^t$ and ${\eta }_{p,f_p}^t$ from the predicted radiation field.

In order to evaluate the effectiveness of proposed MILP model for evacuation problem, several random generated instances are carried out. In addition, the MILP model is coded in YALMIP of MATLAB and solved by Gurobi with a time limit of 600 s. All the scenarios are carried out on a machine with Intel Core i5-7500 3.40 GHz CPU and 8 GB of RAM.

The details of five instances for testing are shown in

Table 2. Evacuation instances.

Instance	Depots	Pickups	Shelters	Buses	Capacity
1	1	5	2	20	25
2	1	5	2	25	20
3	2	10	3	8, 12	25
4	2	10	3	5, 10	25, 30
5	2	10	3	10, 10	25, 30

. All instances ensure that all evacuees can be transferred to shelters by buses. The differences between the first two instances lies in the number and capacity of buses. For complex evacuation network with more nodes, the simulation results are given by the last three instances. Note that the number and capacity of buses might be vary between depots.

The simulation results of five instances are shown in

Table 3. Evacuation simulation results.

Instance	Gap	Compt Time (t)	Evac Time (t)	Radiation
1	0%	22.77	611.48	$1.36\times {10}^7$
2	0%	35.80	611.48	$1.33\times {10}^7$
3	0%	144.98	453.32	$1.63\times {10}^7$
4	0%	231.73	473.42	$1.55\times {10}^7$
5	0.65%	time limit	449.09	$1.53\times {10}^7$

and compared in four categories, optimality gap, computing time, total evacuation time, and total radiation exposure. Note that 4 out of 5 instances reach the optimal solution within the time limit. Instance 5 with large-scale network fails to reach the optimum, and its optimality gap is $0.65\%$. It can be seen from the table that the total evacuation time and the total radiation exposure decrease while the computing time increases with the increase in network size.

The detailed evacuation routes of Bus 7 in Instance 2 and Bus 8 in Instance 3 are illustrated in Figure 6 and Figure 7. The yellow nodes, red nodes, and green nodes represent depots, pickup points and shelters. Obviously, the bus can arrive at different pickup points during different trips to carry out evacuation plan.

Table 4. Evacuation routes of buses in Instance 4.

Bus	Trips				Evac Time
	1	2	3	4
1	8	9	10	4	373.82
2	9	9	10	5	351.21
3	8	10	10	5	355.96
4	8	10	10	5	355.96
5	9	9	10	5	351.37
6	2	2	4	5	473.42
7	3	1	4	5	472.14
8	3				99.22
9	7	6	6	5	465.18
10	7	6	6	5	465.18
11	8	9	8	5	441.09
12	2	1	3	3	449.09
13	7	6	6	4	469.59
14	3				99.22
15	7	6	6	5	465.18

shows the detailed evacuation routes of buses in Instance 4. It can be seen from the table that 13 out of 15 buses run four trips to different pickup points, and there are three buses running the same evacuation routes of {7,6,6,5}. Moreover, the maximum evacuation time of Bus 6 represents for the total evacuation time.

6. Conclusions

In this paper, we have presented a nuclear accident emergency response system based on UAV and bus collaboration. As for the UAV-based radiation field estimation, a Gaussian puff model is adopted to describe the nuclear radiation field. Then, the radiation field estimation is transformed into the parameter estimation problem for diffusion model. Based on the radiation measurement model with white Gaussian noise, a CRLB-based metric is proposed to evaluate the amount of information contained in measurements, which depends on the measurement coordinates. With consideration of the mobility of UAV, the coordinate optimization combined with UAV routing is formulated as an MINLP problem. Through the discretization of time-varying parameters, the MINLP problem is solved by a two-stage solution procedure based on GA at each time interval, and then MLE is applied to estimate the parameters. As for the bus-based emergency evacuation, the radiation dose per second for each trip is assumed to be constant and obtained from the predicted radiation field. Based on the MILP model considering the total evacuation time and the radiation exposure to evacuees, the optimal operating strategies for buses are obtained while ensuring the safety of evacuees. The simulation results demonstrate the effectiveness of our proposed system from both radiation field estimation and evacuation planning. Overall, it can be concluded that the proposed system can accurately estimate the radiation field with limited resources and efficiently carry out nuclear emergency evacuation in a safe manner.

Future extension of this work includes improving the optimization models and developing fast solution algorithm for evacuation planning. For instance, the constraints on the mobility of UAV are extremely strict in coordinates optimization problem, which leads to loss of information contained in measurements. Therefore, the improvement of the constraints on UAV routing is a future extension. Furthermore, since the assumption that the radiation dose is constant for each trip is inaccurate, the evacuation planning in dynamic radiation field is an interesting extension of this work. At last, fast and reliable solution algorithm can also be further researched according to the real-time decision-making requirements of nuclear emergency response.