Mathematical aspects of remote assessment of the radiation state of contaminated areas

The use of radioactive materials is widespread in scientific investigations and various sectors of the economy. There are also extremely radiation-hazardous objects, for instance the well-known Chornobyl Exclusion Zone (Chornobyl, Ukraine) covering the large contaminated areas and the Shelter Object containing the materials of huge radioactivity of about 20 MCi. To safe handling with such objects and materials, the correct their monitoring, detection and characteristics evaluation are vital. The modern development of small flying machines, measurement equipment, and information technologies allow one to increase the amount of measurement data and their accuracy, and to reduce the processing time. On the other hand, the requirements to accuracy, quickness, and correctness of data interpretation increase as well. To solve these problems effectively, the mathematical tools of data processing should be improved. The main mathematical problem at the remote evaluation of radioactive fields relates to the solving the inverse problem for the Fredholm integral of the first kind. In this research, we consider the reconstruction of surface density of gamma radiation on the ground using the data of aerial shooting. We survey the methods for solving the inverse problem, their advantages and disadvantages. The adaptation of the methods to the reconstruction of nonstationary discontinuous radioactive fields is presented. We modify the numerical algorithms using the opportunities of modern calculating software. In particular, it is considered the task when the algorithm reconstructs the density distribution very well.


Introduction
The energy of radioactive decay [1] can be useful supporting the high level of human societies development but conventionally also possesses unpleasant often dangerous effects.
In few recent decades the issues of radiation safety are becoming increasingly important. This is caused by the various reasons [2] including essential increasing the nuclear waste, expansion of the areas of contaminated land, intensive use of nuclear materials in military [3] and civilian purposes, operation of power plants containing the nuclear reactors (nuclear plants, submarines, IOP Publishing doi: 10.1088/1755-1315/1049/1/012015 2 missiles, and etc.), increasing a number of terrorist attacks and level of their preparation. Moreover, due to the prevalence of radiation materials, negative consequences of technogenic and natural disasters (Chornobyl and Fukushima nuclear accidents [4]) can be multiplied.
To organize the effective control and correct handling of radioactive materials and to provide the radiation safety, small flying machines or so-called unmanned aerial vehicles (UAVs) began to apply widely. The advantages [5][6][7][8][9][10][11] of this approach are well known. Among them are: • rapid gathering of data extracted from the large areas; • in contrast to the manned vehicles, UAVs do not require special strict regulations or permissions, there are less problem with pilot's safety; • UAVs can provide more detail explorations of the areas with complicated reliefs or obstacles [2,12,13]; • their exploitation is less expensive; • they are less dangerous for environment (accident does not lead to substantial injury or damage of surrounding structures); • there is a possibility to use them in particularly unfavorable conditions (strong radioactive, magnetic, temperature fields); • in general, UAVs can get as close to the source as possible and, in turn, decrease the influence of attenuation effects; • multipurpose application (relatively easy to change equipment).
When radiation fields are considered, special attention is paid to the development of procedures for identifying the array of detector's readings [8,9,14,15]. The choice and details of the implementation of such procedures depend on the problem's specifics, but, as a rule, they are based on the statistical nature of the signals being recorded. The latter leads to the need to solve the inverse integral problem, which belongs to the class of ill-posed problems and can demonstrate the multi-valudness of solutions or instability. Starting with the works of Tikhonov and Lavrentyev, with the subsequent contribution of many others [16,17], a number of algorithms [18,19] have been proposed to regularize the problems including variational, Friedman, weighted singular schedule.
It is intensively developed the methods for solving the inverse problem arising in statistical statement [2,20,21], i.e. the problems which do not possess exact solutions due to containing, for instance, data with measurement errors or the problems in which the amount of in-coming data is restricted, while a model requires the assessment of huge number of model's parameters.
However, their implementation for a specific task is not easy, in particular in the case of multidimensional, singular or dynamic problems that arise during aerial photography, visualization of areas of high radiation.

Research aim and objectives
It is important the increase of the algorithm's stability, the choice of the optimal procedures for a particular task, ensuring the accuracy of reproduction of the characteristics of the studied system, speed and ease of obtaining results. The aim of this report is to consider integral direct and inverse problems arising in the task of remote monitoring of distributed radioactive fields. In this research we investigate the problem of remote sensing of ground surface to reveal the distributed gamma ray emission and estimate its characteristics.

Statement of the problem
Let us start from the mathematical description of the problem omitting some technical and physical details [9,[12][13][14]. In particular, the simplified scheme of remote sensing is depicted in mounted on a small flying machine providing us the counting of gammas produced on the ground surface. The detector's eyeshot is formed by the collimator with rectangular cross-section and the aperture 2θ rad. Hence, on the ground surface it is appeared the view window (almost rectangular), the gammas from which reach the detector. It is obvious that the size of the window depends on the detector's altitude h d and θ, i.e. the window side 2∆x = 2(h d − h) tan θ, where h is the relief altitude. For instance, when the detector with θ = π/3 rad is located at a height of 55 m above the ground surface, the view window size is 2∆x = 190 m.
Next, the detector traps gammas during the time interval [t i−1 ; t i ]. We can also assume that the exposure time of the detector is t i − t i−1 = τ = const. At the end of this time, the number of gammas W is fixed. Note that, during the time τ , UAV moves with the velocity v m/s along the horizontal axis and the view window is shifted on the distance vτ m as shown in figure 1, i.e. when τ = 1 s and v = 100 km/h, the window displacement is 28 m. As a consequence, if UAV flies the distance about L = 4000 m, we obtain L/(τ v) ≈ 140 values of W . Figure 1. The simplified scheme of remote sensing.
Thus, in general, to evaluate the number of trapped gamma rays at a time moment t i , it can be used the following integral [14,19] equation where w is the density of gammas, x d ∈ R n is the detector's position and x ∈ R n is a point in the domain Ω representing a radioactive material. The quantity where σ is the density of gamma sources; X d is the detector position; s is the effective area of detector's inlet; µ = 1/ℓ is the coefficient of attenuation in air (ℓ stands for the mean free path of the gamma rays [22]); the function g = (h d − h)/r coincides with the cosine of angle between the vertical and the direction from the detector to the point x; r is the distance between detector and the point x lying on the surface. Thus, two problems arise: the former is to estimate the number of events W registered by the detector for a given distribution of σ (direct or forward problem); the latter is to restore the distribution σ, when the data W are given (inverse problem).
In this report, we consider the one-dimensional problem (1). Thus, from equation (1) and taking into account that t i = iτ it follows where Using the substitution vt = x d , relation (2) can be written in equivalent form Thus, below we are going to solve the direct problem evaluating the integral in (2).

Direct problem
To solve the direct problem represented by relation (2)   Note that integral (3) can be calculated by means of conventional rectangle method. Recall that, to approximate the integral of f (x) over the interval (a, b), it can be used the following version of the rectangle method: where ∆ 1 = 2∆x/N 1 , ∆ 2 = τ /N 2 . The function G reads as follows Comparison of specified density σ(x) and results obtained with the in-built function NIntegrate... and via the relation (6) is depicted in figure 2. The profile of W (iτ v) is stretch (each W value is multiplied by 60) over the vertical axis. In this figure, we can see that the doublet structures are distinguished on the resulting profile of W .
Analyzing the structure of the right part of relation (6), it is worth to note that the discretizations over t (adds with k) and x (adds with j) are, to some extent, equivalent. Thus, we can take N 2 not small but arbitrary, whereas N 1 can be chosen 1. Plotting the profile of W at N 2 = 7 and N 1 = 1, we obtain the indistinguishable W profiles.

Inverse problem
The solution of the inverse problem is related to the reconstruction of the function σ, when the values of W (i) are known. As it was mentioned above, such problems are more complicated [17,18] in comparison with the direct problems due to their ill-posedness. To examine these problems, it has been developed a wide range of approaches but there are examples when we can avoid the ill-posedness, as in our case.
In the case of equation (2), we find the function σ(x) in the discrete form. Namely, the rectangle method gives us the natural way to discretize the spatial variable x. Before the application of relation (6), it is worth to note that the argument of σ is not suitable for using since it depends on many parameters. Thus, let us assume that τ v = ∆ 1 N 2 and rewrite it in the following form Since i, N 2 , k, j are natural numbers, then the set of unknown values σ(x) is defined for the multiple values ∆ 1 . Next, putting σ([(i − 1)N 2 + k + j − 1]∆ 1 ) = σ q and now q = (i − 1)N 2 + k + j − 1 are natural, the relations (6) can be presented in the form which is the linear system of algebraic equations with respect to σ q . The function G is defined by (7).
It is obvious that the equations for i ∈ [135; 140] contain the extraquantities σ 141,...,146 . All these unknowns are assumed to be zero. So, using the in-built Mathematica functions, the numerical procedure is constructed. At first, we specify the vector of data W writing them into the array Y =  figure 3, it is worth to notice that there is a point at the right edge lying slightly beyond the theoretical curve σ(x). It is obvious that this does not influence on the other points.

Concluding remarks
Thus, in the research presented above we considered the mathematical issues related to the remote sensing of ground surface. Modern achievements of software, numerical methods, UAVs provided additional opportunities to improve and understand more deep the problem we encounter during remote sensing. In particular, using the tools offered by the software package Mathematica, we constructed the solutions of direct and inverse problems in a simple way. It turned out that the solution of the direct problem conveys very well the localized sources of emission including of doublet type. Furthermore, in these studies we developed the procedure for the inverse problem solving which restored the theoretical distribution perfectly. It is possible that this is related to the specific choice of the quantities h d , h, v, and τ . But this issue requires additional studies.
Presented and similar investigations are useful to consider in connection with other ways of monitoring of the territories affected by accidents at nuclear facilities such as at Chornobyl or Fukushima NPPs. Huge amount of contaminated materials, complex contamination of areas, large temporal scales of radioactive decay highlight novel regularities the nature's evolution and adaptation [1,4,23,24]. A systematic approach to solving existing problems and preventing new ones is the key to sustainable development of human society. The findings presented can be useful for the problems of revealing and specification of weak ionic radiation occurring in homeland security, ecology, nuclear medicine, and etc.