Mathematical tools of solving the problem of restoring the surface distribution of radiation pollution based on remote measurement data

The modern achievements in the construction of small flying machines cause the active development of remote monitoring, in particular geophysical airborne gamma-ray spectrometer surveying. Such observations are important, since large amounts of man-made radioactive materials get into the environment, especially during accidents like at the Chornobyl or Fukushima nuclear plants. On the other hand, the natural distribution of radioactive sources is inhomogeneous and can provide us with useful information about the soil structure. One of the problems appearing at the handle of information collected with unmanned aerial vehicles concerns the correction of readings to identify the peculiarities of gamma-ray fields. To perform this, the analytical method based on the solution of the inverse problem formulated in terms of integral relation is used. In this research, to reconstruct the surface distribution of the gamma-ray field, the Tikhonov and Landweber techniques are applied. It is shown that these algorithms allow one to distinguish radioactive hot-spots located closely.


Introduction
To support the sustainable development of economies and human societies, the effective implementation of innovative technologies, methods of rational nature management, and safety components should be realized.To a large extent, this applies to handling the radioactive materials which can provide significant benefits and at the same time can cause substantial damage.To organize the effective control of radioactive materials accompanied by ensuring radiation safety, the novel achievements in the fields of remote sensing and information technologies [1][2][3][4] can be useful.There are many different means for gathering and processing information on observed objects including satellites, manned aviation, ground observation, and 1254 (2023) 012099 IOP Publishing doi:10.1088/1755-1315/1254/1/012099 2 etc. [5].Among them, the small flying machines or so-called unmanned aerial vehicles (UAVs) should be noted, the use of which has developed rapidly due to their advantages [1,[6][7][8][9] including the possibility to reveal the size-limited gamma-ray anomalies [10].
Nevertheless, the collected data undergo the influence of different effects, i.e.UAV velocity and high, detector's properties, attenuation, natural radiation, and require further processing and corrections.To improve the quality of UAV measurements, it is utilized the approach [11][12][13][14][15] which is based on the use of a solution of inverse problem.This problem involves a multiple integral relation, approximation of which, in turn, causes well-known problems.However, the main obstacle relates to the construction of a solution to the inverse problem, since we deal with an ill-posed problem that produces strong instability, oversmoothing effects, and etc.To overcome this difficulty, many techniques have been developed including the Tikhonov and Landweber algorithms, which we used in these studies.

Research aim and objectives
In this research, we consider the adaptation of algorithms for the inverse problem which arises during the remote monitoring of distributed radioactive (namely gamma-ray) fields.In particular, restoration of the gamma field, distributed over the ground surface, and identification of its peculiarities (local intensities) are performed.

Statement of the problem
As known from numerous experiments [15][16][17], the readings of airborne detector require some corrections.One of the possible ways to do this is the inversion of airborne data using the algorithm of the proper inverse problem.The mathematical statement of this problem involves the integral relation describing the attenuation of the gamma radiation with distance from the ground surface, the distribution of radioelement sources in the ground, and the detector's properties.The UAV velocity also should be taken into account.
Thus, consider the simplified scheme (figure 1(a)) of remote sensing, continuing our studies presented in [6].A small flying machine is equipped with a detector providing us the counting of gammas emitted from the ground surface.UAV moves along the axis Ox with the velocity v.It is marked by the point D with the coordinates (x d , y d ) and the altitude h d in figure 1(a).
Let us assume that the detector's eyeshot is formed by the collimator with a rectangular cross-section and the aperture 2θ rad.Then the view window appearing on the ground surface is almost rectangular.The size of this window depends on the detector's altitude h d and θ.The elementary geometric reasoning leads us to the simple relation for the window side 2∆x = 2∆y = 2(h d − h) tan θ, where h(x, y) stands for the relief altitude.Dealing, for instance, with the detector with θ = π/3 rad located at a height of 55 m above the ground surface, the view window size is 2∆x = 190 m.
During a flight, the detector collects gammas and fixes the number of gamma quants W (t i ) at discrete moments of time t i , i = 1, 2, . . ., N .For the sake of simplicity we will assume that all time intervals [t i−1 ; t i ] are of the same length and, therefore, the detector's exposure time is t i − t i−1 = τ = const.Thus, if we fixed τ = 1 s and v = 27 m/s (about 100 km/h), then UAV flying L = 4000 m provides us with L/(τ v) ≈ 140 ≡ N values of W .
Thus, in general, to evaluate the number of trapped gamma rays at a time moment t i , it can be used the following integral [3, 11-13, 15, 16] where σ is the density of gamma sources; X = (x, y) ∈ R 2 is a point in the domain Ω representing a radioactive material; Λ is the scaling factor depending on the geometric and physical properties of detector.The quantity where µ is the coefficient of attenuation in air; the function φ = (h d − h)/r coincides with the cosine of angle between the vertical and the direction from the detector to the point X is the distance between detector and the point X lying on the surface.Thus, the problem arises: to restore the surface distribution σ(x, y), when the data W are given.This is a well-known inverse problem.

Construction of the algorithms for the inverse problem
Before solving the inverse problem, we prepare the test data W for checking the quality of the algorithm's work.In fact, evaluating the quantities W (t i ), the forward problem is solved.Thus, we specify the function of radioactive sources density where U nitBox[z] stands for the in-built Mathematica package function, equal to 1 for |z| ≤ 1/2 and 0 otherwise.Using the function U nitBox[z], we model the localized sources of high radioactivity.In other words, four peaks on the density function σ can represent, for instance, trenches containing radioactive waste.Actually, in this research, the function (2) doesn't depend on y and its graph is depicted in figure 1(b).The model (1) also requires the parameter µ = 0.001 1/m; relief profile h(x, y) = 20 m; UAV altitude h = 75 m.Although we can incorporate in relation (1) the UAV specified velocity v i at each interval [t i−1 , t i ], but in this research we neglect small deviations of UAV velocity and it is assumed that the UAV moves all flight time with the mean constant velocity v = const.Since the UAV's movement is performed along the axis Ox, its coordinate y d = 0.
Then, substituting the new variable u = vt and other parameters and functions into relation (1), we get where For convenience, we take Λ = σ(0, 0)/W 1 , in order to normalize data and to compare with the parent distribution σ.
For instance, if we fix δ x = 9 m, then q = 3, ε = 10, τ = 1 s, v = 27 m/s, N = 144, then structure of the system (6) can be understood from the first two equations We see that the system is the set of coupled equations.The matrix K is multi-diagonal and not square, in general.This means that the system does not possess the unique solution and as a rule it is ill-conditioning.To obtain the proper solution in this case, the regularized algorithms are used.In this research, we apply the Tikhonov and Landweber algorithms [18].

Tikhonov method
Let us recall shortly that the Tikhonov regularization technique is based on the singular value decomposition (SVD) procedure applied to the matrix K and elimination of singular elements from the system's solution which produce the solution's instability.As shown in [18], using the SVD presentation of the m × n real-valued matrix K in the form K = U ΣV , where U and V are the m×m and n×n unitary matrices respectively, Σ = diag(λ i ) is the m×n diagonal matrix with diagonal nonnegative entries λ i called singular values of K. Then the solution of system (7) can be cast in the matrix form To get rid of small singular values of the matrix diag(λ −1 i ), the Tikhonov filter function transforming the small singular values is used.Thus, instead of exact system's solution we have Using the identity , where I is the unit matrix, the preceding relation reads as follows where α is a regularization parameter.Thus, at first consider the properties of matrix K used in this research.Instead of direct evaluation of K's singular values, it is simpler to evaluate eigenvalues of K K using the in-built Mathematica function Eigenvalues[].The resulting distribution of eigenvalues λ 2 i is presented in figure 2(a) (here the semilogarithmic plot is used) and shows the presence of small singular values for large i.Finally, putting α = 10 −3 , the solution ( 8) is evaluated and after rearranging ((−20+2j)hx/2; (g α ) j ), j = 1, . . ., N the figure 3(a) is plotted.When α decreases up to 0.5•10 −8 , the solution ( 8) is depicted in figure 3(b) and the growth of spurious fluctuations is observed.It means that α is too small.
In these studies we also are interested in the solution when the components of system (7) are desturbed by noise modeled with the random normally distributed numbers ξ.Thus, we consider the system (7) in the form w = (K +ξ)g, where ξ ∼ N orm(0, b).If we take the standard deviation b = 5 • 10 −4 , the Mathematica function RandomVariate[NormalDistribution[0, b], Nq] generates N q random numbers.The random data histogram is depicted in figure 4(b) and tells us about a good approximation of normal distribution.As expected, the corresponding solution of inverse problem plotted in figure 4(a) shows more irregularity in comparison with the non-noised solution from figure 3. When b = 0.001 and random numbers undergo a greater deviations from their mean value (figure 4(d)), the solution (8) drawn in figure 4(c) becomes even greater irregularity, so that doublet structures become faintly visible against the background of fluctuations.The similar solution's behaviour is observed, when noise is added to the vector w in system (7) but at the smaller values of b.

Landweber method
This technique belongs to the iterative regularization methods and is written in the form of iterative procedure [18] where , η is a scalar controlling the speed of iteration convergence.As shown in [18], the procedure (9) can be reduced to the matrix representation Let us first fix not large iteration number ν = 50.Then the corresponding solution (10) is shown in figure 5(a) and represents a good discrimination of the doublet pairs.
When the number of iterations are sufficiently large, i.e. ν = 500, the solution ( 9) is plotted in figure 5(b).We see that the difference between exact solution and its approximation becomes smaller.It should be noted that the development spurious oscillations at the left edge of the interval, whereas the solution at the right edge behaves well.The reason for this instability requires additional studies, although this effect does not influence the hot-spot detection.

Concluding remarks
Thus, the research outlined above dealt with the mathematical issues devoted to the remote sensing areas contaminated by the highly radioactive materials.Nowadays the use of UAVs for remote sensing grows rapidly, therefore, support and assistance of UAV technology are extremely important.To do this, modern software and numerical methods are utilized.In particular, the software Mathematica package offers the tools for constructing the solutions of forward and inverse integral problems in a simple way.Using this package, we adopted the Tikhonov and Landeweber algorithms for the inverse problem and applied them to the test problem which takes into account the surface density distribution of gamma sources with doublet structures, while other quantities were chosen as simply as possible.It turned out that the algorithms for the inverse problem reconstruct the close doublet structures quite well, remaining them separated.
Presented studies can be used for organizing the monitoring of the territories with man-made radioactive waste, for instance affected by accidents at nuclear NPPs [14].As is well known [6], the territory of the Chornobyl Exclusion Zone contains huge amounts of spent nuclear fuel, radioactive waste, and its temporary localization.These storage facilities consist of complex "Vector", radioactive waste disposal points "Buryakivka", "Pidlisnyi", "ChNPP Stage III", nine temporary radioactive waste confinement sites with a total area about 10 hectares and include nearsurface disposal facility, trenches, and mound storages.In order to distinguish the contours of lost trenches, the use of presented findings is promising.
The results obtained can be interesting for the experts in the fields of remote sensing of radiation, homeland security, ecology, nuclear medicine, and etc.

Figure 1 .
Figure 1.The simplified scheme of remote sensing (a) and the surface of the function σ(x, y)the density of gamma sources (b).
Figure 2(a) represents the comparison of the exact graph of the density function σ and data W .

Figure 2 .
Figure 2. The comparison of the profile of the function (2) and data W evaluated by the formula (3). b