%\setcounter{chapter}{3}
\chapter{Simulation of proposed algorithms}

In order to find parameters that influence applicability of proposed algorithms, a~series of simulations are performed. The simulations are not complete – among other things, kinematics of the robot is not considered because it is not important in relation to the localization algorithms. Responses of detectors are simulated only in points of interest, for example a~movement between lines is not considered in case of mapping algorithm. System MATLAB has been used to realize simulations.

\section{Simulation of detector response}

A model of a~source of ionizing radiation may be described by vector $\bm{\theta}_r =$ \vect{x_r, y_r, \lambda_r} where $x_r$ and $y_r$ are coordinates of the source and $\lambda_r$ is its intensity represented by mean emission of particles (photons in case of gamma radiation). Given $R$ sources, individual models can be joined into a~single vector $\bm{\theta} = [\bm{\theta}_1, \bm{\theta}_2, \cdots, \bm{\theta}_r]^\text{T}$. Radioactive decay of a~source is a~random process that can be described by Poisson distribution \cite{poisson}. The probability of emission of $x$ particles is given by:

\begin{equation}
p(x=X) = \mathcal P (x; \lambda) = \frac{e^{-\lambda}\lambda^x}{x!}
\end{equation}

Making a~simplifying assumption that conversion efficiency of detectors is 100~\%, the number of emitted particles correspond to measured total count. 

Another source of particles that has to be considered is natural radiation background. Since the background comes not only from natural radionuclides but also from cosmic radiation and other sources, it is not feasible to model it precisely. It has been empirically discovered that contribution of radiation background in the area where all experiments took place ranges between 100 and 160 counts per second. For modeling, a~continuous uniform distribution has been used. This distribution is used by function \texttt{rand} in MATLAB. Total count detected by a~detector $k$ can be obtained from equation \cite{bayes}:

\begin{equation}
\lambda_k(\bm{\theta}) = \lambda_B + \sum_{r=1}^R c_k(\bm{\theta}_r)
\end{equation}
where $\lambda_B \leftarrow \mathcal U([100, 160])$ is the contribution by background, and $c_k(\bm{\theta}_r)$ is a~count rate due to source $r$. Its theoretical value is given by integral \cite{bayes}:

\begin{equation}
c_k(\bm{\theta}_r) = a_r \int_0^T\frac{\text{d}t}{||\bm{x}_k - \bm{x}_r + \bm{\dot{x}}_kt||^2}
\end{equation}
where $a_r \leftarrow \mathcal P(\lambda_r)$, $T$ is sampling period, $\bm{x}_k$ is position of detector $k$, $\bm{x}_r$ is position of source $r$, and $\bm{\dot{x}_k}$ is velocity of detector $k$. Since the $\lambda$ values are typically in order of thousands and the Poisson distribution is numerically stable for order of tens at most, it is necessary to use approximation by Normal distribution \cite{normal}. Therefore $a_r \leftarrow \mathcal N(\lambda_r, \lambda_r)$. The integral may be solved analytically but it provides a~model which is too much precise. Sufficient approximation is provided by formula:

\begin{equation}
c_k(\bm{\theta}_r) = \frac{a_r}{||\bm{x}_k(t+T) - \bm{x}_r||^2}
\end{equation}

Another effect that has not been discussed yet is a~directional characteristics of the detector. That has been measured and presented in Sec.~\ref{sec.dirsens}. Height of the detector above terrain $h$ should be considered as well. The previous equation therefore has to be updated to:

\begin{equation}
c_k(\bm{\theta}_r) = \frac{K_k(\phi_{k,r}) a_r}{||\bm{x}_k(t+T) - \bm{x}_r||^2 + h^2}
\end{equation}
where $K_k(\phi)$ is sensitivity in direction $\phi$, and $\phi_{k,r}$ is an angular coordinate of source $r$ in a~coordinate system of detector $k$. That can be calculated as:

\begin{equation}
\phi_{k,r}=
  \begin{cases}
    \arctan\frac{y_r'}{x_r'},& \text{if } x_r' > 0\\
    \arctan\frac{y_r'}{x_r'} +180^{\circ},              & \text{otherwise}
  \end{cases}
\end{equation}
where

\begin{subequations}
\begin{align}
x_r' = (x_r-x_k)\sin{\alpha} + (y_r-y_k)\cos{\alpha}\\
y_r' = (y_r-y_k)\sin{\alpha} - (x_r-x_k)\cos{\alpha}
\end{align}
\end{subequations}
where $\alpha$ is an azimuth of the robot carrying the detector. Assuming constant velocity between two sampling points, coordinates of detector  $\bm{x}_k = [x_k, y_k]$ in time $(t+T)$ are given by:

\begin{subequations}
\begin{align}
x_k(t+T)&= x_k(t) + ||\bm{\dot{x}}_k(t)||T \cos(90^{\circ} - \alpha)\\
y_k(t+T)&= y_k(t) + ||\bm{\dot{x}}_k(t)||T \sin(90^{\circ} - \alpha)
\end{align}
\end{subequations}

\section{Mapping algorithm}

Mapping of the area of interest is simulated first. For both simulations and real experiments, a~site near Faculty of Electrical Engineering and Communication in Brno has been chosen. An ortophoto of the site with the marked polygonal area is shown in Fig.~\ref{fig.ortophoto}. A position of the base station is displayed there as well. The chosen polygon has area of 438 m$^2$.

\begin{figure}[!h]
  \begin{center}
    \includegraphics[width=12cm]{images/qgis-orto.jpg}
  \end{center}
  \caption{Ortophoto of the experimental site \cite{holomas}}
  \label{fig.ortophoto}
\end{figure}

A result of simulation using line separation of 1 meter can be seen in Fig.~\ref{fig.simmap0}. Function \texttt{scatter} is used to display individual measured data points. Theirs color corresponds to a~sum of counts per second measured by both detectors according to color bar displayed in the figure on the right. Logarithmic scale is used because CPS values are spanned over multiple decades. An interpolated map is presented using fuction \texttt{contourf} in Fig.~\ref{fig.simmap1}. It is obvious that one strong and two weak sources are present in the polygon. Theirs estimated positions are displayed by using black crosses.

\begin{figure}[!h]
  \begin{center}
    \includegraphics[height=9.9cm]{images/sim_mapping0.pdf}
  \end{center}
  \caption{Simulation result of mapping}
  \label{fig.simmap0}
\end{figure}

\begin{figure}[!h]
  \begin{center}
    \includegraphics[height=9.9cm]{images/sim_mapping1.pdf}
  \end{center}
  \caption{Interpolated map}
  \label{fig.simmap1}
\end{figure}

In order to find parameters for estimation of the position, the simulation was repeated fifty times for each value of parameter within defined range. An optimal value for parameter can be selected using conveniently chosen criterion. First, parameter $k$ mentioned in Sec.~\ref{sec.interp} was sought. The following criterion is defined to evaluate optimality:

\begin{equation}
J_1 = |n_E - n|
\end{equation}
where $n_E$ is estimated count of sources, and $n$ is the actual count. Dependence of criterion $J_1$ on value of parameter k is shown in Fig.~\ref{fig.map_crit} (left). The optimal value is 0.996. Second parameter to be found is a~radius $r$ within which data points should be interpolated in order to get more precise estimation of the position. Sum of distances between actual and estimated positions is used as optimality criterion:

\begin{equation}
J_2 = \sum_{r=1}^R \sqrt{(x_{E,r}-x_r)^2 + (y_{E,r}-y_r)^2}
\end{equation}

Dependence of criterion $J_2$ on the value of parameter $r$ can be seen in Fig.~\ref{fig.map_crit} (left). There are two valleys  in the dependence. The first one probably corresponds with the weak sources, the other one with the strong source. Value 1.5 is chosen as optimal.

\begin{figure}[!h]
  \begin{center}
    \includegraphics[height=6.7cm]{images/mapping_crit.pdf}
  \end{center}
  \caption{Optimality criteria for mapping parameters}
  \label{fig.map_crit}
\end{figure}

\section{Strong source search algorithm}

A possible result of simulation of strong source search algorithm is illustrated in Fig.~\ref{fig.simloc}. Peaks detected in the individual lines are marked by using red crosses. It can be seen in the figure that no peak has been found in the first line, but the next two lines have corresponding peak. That results in alternation of the original trajectory. Part of final loop is displayed as well, and it shows where it is crossing the perpendicular line leading to the source (black cross).

\begin{figure}[!h]
  \begin{center}
    \includegraphics[height=9.9cm]{images/sim_locating.pdf}
  \end{center}
  \caption{Simulation result of the strong source search}
  \label{fig.simloc}
\end{figure}

Parameter $d$ needs to be set. It defines maximal allowed distance of neighboring peaks in which the trajectory is altered. Optimality criterion is defined as a~multiple of the number of basic lines $l$ and distance between position of the source and the loop crossing.

\begin{equation}
J_3 = l \sqrt{(x_S-x_L)^2 + (y_S-y_L)^2}
\end{equation}

Dependence of criterion $J_3$ on parameter d is shown in Fig.~\ref{fig.loc_crit}. The optimal value is ambiguous, eventually 1.4 was chosen.

\begin{figure}[!h]
  \begin{center}
    \includegraphics[height=6.7cm]{images/locating_crit.pdf}
  \end{center}
  \caption{Optimality criterion for a~peak distance parameter}
  \label{fig.loc_crit}
\end{figure}

\begin{figure}[!h]
  \begin{center}
    \includegraphics[height=9.9cm]{images/sim_circle1.pdf}
  \end{center}
  \caption{Simulation result of the first circular algorithm}
  \label{fig.simcirc1}
\end{figure}

\section{First circular algorithm}

A result of the simulation of the first circular algorithm is shown in Fig.~\ref{fig.simcirc1}. The circle has radius of 4 meters. After the data along circle are measured, direction of the source may be estimated. Accuracy of estimation depends on a~choice of parameter $p$ – shift between interpolated and actual maximum. As optimality criterion a~distance of the source and a~line passing through circle center and directional point  is used.

\begin{equation}
J_4 = \frac{ax_S + by_S + c}{\sqrt{a^2 + b^2}}
\end{equation}
where \vect{a,b,c} are parameters of the line. Dependence of criterion $J_4$ on parameter $p$ can be seen in Fig.~\ref{fig.circ1_crit1}. The optimal value is 0.7. Once the direction is found, algorithm proceeds to navigation to the source. There are two more parameters that need to be set. One of them is coefficient $K$ which determines how much data from both detectors must differ in order to enable an update of the direction. The second parameter $\Delta\phi$ determines how much the direction is altered. Both parameters have common optimality criterion – distance of loop crossing and the actual source position.

\begin{equation}
J_5 = J_6 = \sqrt{(x_S-x_L)^2 + (y_S-y_L)^2}
\end{equation}

Dependence of criteria $J_5$ and $J_6$ on respective parameters is shown in Fig.~\ref{fig.circ1_crit2}. The optimal value of $K$ is 1.04, optimal value of $\Delta\phi$ is approximately 18°.

\begin{figure}[!h]
  \begin{center}
    \includegraphics[height=6.7cm]{images/circle1_crit1.pdf}
  \end{center}
  \caption{Optimality criterion for a~peak shift parameter}
  \label{fig.circ1_crit1}
\end{figure}

\begin{figure}[!h]
  \begin{center}
    \includegraphics[height=6.7cm]{images/circle1_crit2.pdf}
  \end{center}
  \caption{Optimality criteria for navigation to source parameters}
  \label{fig.circ1_crit2}
\end{figure}

\section{Second circular algorithm}

A difference between this and previous algorithms is only in number of sources detectable from the circle. Data simulated along a~circle compose a~signal shown in Fig.~\ref{fig.circ2_crit} (left). Since the signal is very noisy it may not be apparent that two prominent peaks should be there. Relative peak prominence $p$ which helps to filter those peaks corresponding only to the sources needs to be found. Optimality criterion is simple – a~difference between the expected and detected number of peaks $n_D$.

\begin{equation}
J_7 = |n_D - n|
\end{equation}

Dependence of criterion $J_7$ on peak prominence $p$ can be seen in Fig.~\ref{fig.circ2_crit} (right). The optimal value is 0.8. A complete result of simulation of this algorithms is then shown in Fig.~\ref{fig.simcirc2}. Both sources have been localized successfully.

\begin{figure}[!h]
  \begin{center}
    \includegraphics[height=6.7cm]{images/circle2_crit.pdf}
  \end{center}
  \caption{Finding peaks on a~circle}
  \label{fig.circ2_crit}
\end{figure}

\begin{figure}[!h]
  \begin{center}
    \includegraphics[height=9.9cm]{images/sim_circle2.pdf}
  \end{center}
  \caption{Simulation result of the second circular algorithm}
  \label{fig.simcirc2}
\end{figure}

\section{Third circular algorithm}

This time, the whole area of polygon has to be covered using multiple circles. A possible solution is illustrated in Fig.~\ref{fig.polycirc3}. It is assumed that the sources are detectable 4 meters away from circles. Circle $c_1$ has radius of 3.5 m, other two circles have radius of 4 m. Therefore, the inner area of all circles is covered. Covering of outer areas is sketched using dashed circles. Some small regions are outside the detectable area, but in real conditions the range has not sharp borders. Sources within these peripheral regions would be probably detected as well.

A result of simulation is shown in Fig.~\ref{fig.simcirc3}. It is worth noticing that one source lies very close to circle $c_3$, and for that reason the final loop is omitted – the circle provides additional data for interpolation.

\begin{figure}[!h]
  \begin{center}
    \includegraphics[height=9.9cm]{images/polygon_circ3.pdf}
  \end{center}
  \caption{Layout of circles for the third circular algorithm}
  \label{fig.polycirc3}
\end{figure}

\begin{figure}[!h]
  \begin{center}
    \includegraphics[height=9.9cm]{images/sim_circle3.pdf}
  \end{center}
  \caption{Simulation result of the third circular algorithm}
  \label{fig.simcirc3}
\end{figure}

\section{Fourth circular algorithm}

The polygon may be completely covered by using a~higher number of smaller circles. A possible solution is illustrated in Fig.~\ref{fig.polycirc4}. As in the previous section, it is assumed that sources can be detected within 4 meter radius. Three larger circles have radius of 3 m, the others have radius of 2 m .

First, only circle data are simulated in order to test algorithm for filtering intersections of directional rays. A result of this simulation with marked rays is shown in Fig.~\ref{fig.simcirc4rays}. There are two intersections inside circle $c_5$. The intersection that is nearer to the source is accepted, the other one is discarded according to presented rules. Three intersections that are near to the source in the middle of the figure are united into a~single point. Choice of respective parameter $D_{min}$ is not critical here because all sources are well spatially separated. Its value may be set for example to 1.5 m. The last visible intersection is caused by false positive detection of the peak on circle $c_2$. Because the corresponding intersection is discarded, the ray would be classified as a~potential source with known direction and the robot would have tried to find it inside circle $c_2$. The last source have only one matching ray, and is therefore also classified as a~potential source with a~known direction.

A complete result of simulation can be seen in Fig.~\ref{fig.simcirc4}. The black crosses represent the actual position of the sources, the red crosses and the red ray represent estimations. In order to maintain the figure well arranged, search of detected non-existing potential source is not displayed -- it would have added more data points inside respective circle. The upper left potential source is being approached from the circle, the other two sources are being approached from a~point where localization of a~previous source have ended. Only a~part of data points is visualized -- those in which update of direction is allowed.

\begin{figure}[!h]
  \begin{center}
    \includegraphics[height=9.9cm]{images/polygon_circ4.pdf}
  \end{center}
  \caption{Layout of circles for the fourth circular algorithm}
  \label{fig.polycirc4}
\end{figure}

\begin{figure}[!h]
  \begin{center}
    \includegraphics[height=9.9cm]{images/sim_circle4_rays.pdf}
  \end{center}
  \caption{Simulation of direction rays}
  \label{fig.simcirc4rays}
\end{figure}

\begin{figure}[!h]
  \begin{center}
    \includegraphics[height=9.9cm]{images/sim_circle4.pdf}
  \end{center}
  \caption{Simulation result of the fourth circular algorithm}
  \label{fig.simcirc4}
\end{figure}