Method of Multi-Angle Transmission Radiowave Tomography of Dielectric Objects

: A method for solving the inverse problem for reconstructing the spatial distribution of dielectric permittivity from the results of multi-angle transmission broadband radiosounding is proposed. The method is based on inverse wave propagation. The average refractive index of the medium along the wave trajectory is calculated by comparing the results of the calculation of the time delay of the inverse signal in the entire sounding region and the forward propagation time in a homogeneous medium. This method takes into account diffraction effects in solving a direct problem, which allows one to obtain a resolution in the order of a wavelength. The combination of time delays obtained at different probing angles allows the restoration of the distribution of the refractive index in the medium. The paper presents the results of the numerical simulation of this method. The novelty of the proposed approach compared to the conventional back-projection algorithm is that ray approximation is not applied. Instead of the absorption coefficient (used in X-ray tomography), a time delay is considered, which is restored in the entire probed region. The developed method can be widely used in radiowave tomography or microwave tomography for remote non-destructive testing, diagnostics for the internal structures of inhomogeneous media and the restoration of the shapes of opaque objects based on multisensor sensing. methodology, validation, analysis, investigation, resources, curation, writing—original draft preparation, D.S.; writing— review and editing, K.Z.; visualization, D.S.; supervision, K.Z.; administration,


Introduction
Radiowave tomography [1][2][3][4][5][6] or microwave tomography [7][8][9][10][11][12][13] is used for remote nondestructive testing, diagnostics for the internal structures of heterogeneous media and the restoration of the forms of opaque objects on the basis of multi-sensory sounding. Radiowave tomography differs from X-ray tomography due to the significant impact of wave effects. Since the wavelength of radio waves is comparable with the size of the heterogeneities, wave effects appear, which do not allow the use of the known algorithms of X-ray tomography [14][15][16][17][18][19][20][21], working in the ray approximation. The back-projection algorithm is most commonly used for 3D tomography image reconstruction. Wave effects are possible even for X-ray radiation when applied at the nanoscale, but this problem is solved by switching to higher-energy radiation or charged particles [22]. However, there are approaches to extract phase information to take into account wave effects in the projection approximation [23][24][25][26][27].
Radiowave tomography began development with the availability of digital technologies for the measuring and processing of radio wave signals [1]. In study [1], a tomographic method was proposed for visualizing the reflection and transmission coefficients of the objects under study. To obtain information on the dielectric permittivity, more complex processing algorithms and multiangle measurements are required.
Transmission radiowave tomography involves irradiating the object under study on one side and measuring the spatial amplitude-phase distribution behind the object. The source and the receiver are located on opposite sides of the probing object. In this method, the recorded field is blurred and has no sharp boundaries, as in the case of X-ray tomography [14,15], because of the diffraction effect. The effects of refraction and multiple reflection remain even if the ray approximation is appropriate. When developing methods of tomography using radio waves, it is necessary to take into account the wave nature of such radiation and the phase delays of the signal in the medium. Taking into account the wave effects requires the development of approaches for solving the direct and inverse problems.
There are different approaches to solving the problems of transmission radiowave tomography. Cross borehole radar methods used to study the medium between the boreholes are widely used [28][29][30][31][32]. In fact, cross borehole radar methods implement transmission radiowave tomography. Approaches are used based on the analysis of the time delay of the signal transmission. The travel time tomography method allows the visualization of dielectric inhomogeneities [28]. However, methods of this kind do not take into account the wave nature of radiation, but ray approximation is applied. Similar approaches are used in seismic tomography [33,34], where the seismic wave velocity distribution in the medium is reconstructed by time delays. The use of such methods is widespread in the method of ultrasonic nondestructive testing [35], in which, similarly to electromagnetic waves, in the scalar approximation, ultrasonic waves are described by the wave equation. However, if the effects of diffraction and refraction are significant in a medium, travel time tomography methods are not effective [36].
Another problem is the resolution of the resulting images in radiowave tomography. To increase the resolution of tomographic images, it is necessary to localize the region of the interaction of radiation with matter. There is already a method of double focusing-the focusing of emitted and received waves with the help of dielectric lenses [37]. As a result, it is possible to form a relatively localized working region, within which the radiation has the form of a collimated wave beam. At the same time, multiple interactions sharply decrease, but the diffraction effects remain significant. Similar principles are used in optical tomography [38,39], when lenses are used to focus optical waves.
In most cases, the entire tomographic system can be considered as linear. Therefore, it can be described in terms of systems of linear algebraic equations, where the parameters of the medium are unknown quantities and the known values are the set of measured signals. To solve such systems, regularization methods are used [31,32,40], since, as a rule, the matrix describing the system has a high conditionality.
There are iterative approaches to solving the inverse tomography problems, which are based on minimizing the difference between the numerical solution of the direct problem and the actual measured data [13]. In [13], a variant of the method of gradient descent is proposed, which makes it possible to reconstruct flat distributions of electrical conductivity. The disadvantages of this method are the computationally complex iterative process and the probability of finding not a global but a local minimum. To avoid falling into the local minimum stochastic, Monte Carlo methods are used [25]; however, there remains the computationally costly necessity of the iterative search of solutions. Nevertheless, in view of the constant development of computer technology, the problem of the large amount of computation becomes less significant.
An effective description of the residual diffraction effects is given by the phase approximation of the Huygens-Kirchhoff method [41], in which the main interaction of radiation with translucent objects is described by a phase screen. The phase distribution in the plane of the phase screen at each point is given in the geometric optics approximation, as if the corresponding wave were flat and propagated strictly along the line through the probing object.
The best results can be achieved by methods based on taking into account the wave nature of radiation [42][43][44][45][46][47] and considering phase information. The use of broadband signals or measurement at multiple frequencies makes it possible to minimize the artifacts of reconstructed tomographic images [42,43] even with unfilled spatial measurements. The radiowave tomographic problem can be solved even in the case of the absence of phase measurements [48] by the application of frequency scanning and the synthetic aperture technique. In [45], the solution of the inverse problem reduces to the method of least squares, but such an approach is quite computationally consuming; however, it potentially allows the achievement of the most accurate solution. For restoring the properties of the medium, the most accurate results allow the obtaining of iterative methods. The optimized search schemes [47] allow the obtaining of a solution within an acceptable time frame. However, there is a need to speed up the solution by searching for the most accurate first approximation for iterative methods.
In this paper, we propose a method of the inverse problem numerical solution for reconstructing the distribution of the refractive index in a medium from the data of multi-angle transmission broadband radio sounding. The proposed approach is new and has not been previously considered. In contrast to the classical method of back-projections with a ray approximation, our method takes into account the propagation of waves of finite length, and the shadow delay of the propagation of waves caused by the presence of dielectric inhomogeneities is considered as shadow projections. The existing methods of diffraction tomography consider cases of slightly contrasting inhomogeneities with an insignificant phase incursion in the inhomogeneities relative to the unperturbed case [49][50][51][52][53][54].
In our case, the contrast of the refractive index can be significant; the ambiguity of the phase change is eliminated by applying measurements in a wide frequency band.
In contrast to holographic tomography, in proposed article, the field diffraction at a variety of frequencies is considered [55]. In this case, the synchronization of the time counting for all frequencies is necessary. We considered multi-frequency sounding. The fundamental difference from existing methods is the use of a time delay as a parameter to determine the refractive index.
Below, we have provided a detailed description of this method. Analytical calculations are given together with the results of the numerical modeling for the clarity and convenience of perception of the material.

Method of Forward Task Solution and Numerical Modeling
It is proposed to consider the scheme of transmission multi-angle sounding by means of radio waves ( Figure 1). The investigated dielectric objects are located on a rotating platform. A broadband plane wave probes objects at each angle of rotation. The amplitude-phase distribution of the field in the plane behind the objects across the set of frequencies is measured. The receiving array is paced at a distance, R, from the center of rotation of the objects under study. Elements in the array are located in steps of less than half the shortest wavelength, that is, the grid is filled. A direct problem is proposed to be solved using the algorithm described in [4]. The algorithm is based on splitting the medium into thin layers and calculating the propagation field in the plane wave spectrum for homogeneous layers of different materials ( Figure 2). Then, solutions for different homogeneous layers would be combined to get the field distribution in heterogeneous media. Similar calculations were made for various angles of rotation of the test objects. As a result, a field distribution was obtained for the given two rectangular objects (Figure 3) in the measurement area, as shown in Figure 4. We modeled additive noise with a uniform statistical distribution in the in-phase and quadrature components of the measured signal at all frequencies to further verify the stability of the method for solving the inverse problem to noise. The in-phase component image (Figure 4a) with the addition of uniform noise with a spread of 30% of the maximum signal amplitude is shown in Figure 4b, and that with the addition of noise with 100% spread, in Figure 4c.
We denote the result of forward task solution by the function: where = the angle of the test object rotation in the XOZ plane, = the cyclic frequency of the probing radiation, and = the position of the measurement system in the z axis.
On the basis of the obtained field in the measurement region, ( , , , ), in the frequency band (in the numerical model from 16 to 20 GHz, 16 frequencies) for different rotation angles, it is necessary to find the distribution of dielectric heterogeneities.
For real experiment function, ( , , , ) needs to be represented as an array of complex numbers that corresponds to the magnitude and phase of the radiowave field passed though the test object and measured at certain points with a fixed spatial step.

Method of Inverse Task Solution
The solution of the inverse problem is based on the inverse propagation of waves from the measurement region towards the emitter. We calculate the field in the sounding region on the basis of the formula: If dispersion effects are not taken into account, then at all frequencies, the delay in matter occurs identically, so after returning to the time domain function, ( , , , , ), it will have the form of a short pulse with a maximum magnitude at = Δ ( , , , ). Thus, finding the maximum magnitude of ( , , , , ) in time, we can find the value of Δ ( , , , ), which describes the time delay of the pulse in a heterogeneous medium relative to that in the homogeneous medium in the entire space. Figure 4 shows the result of calculating the delay function, Δ ( , 0, , ), for the test objects shown in Figure 3a. In Figure 5, darker areas correspond to larger values of delay. Additionally, in Figure 5, it can be seen that an object with a large refractive index corresponds to a darker color in the figure, or larger delay values. This delay is proportional to the product of the refractive index minus 1 and a distance (2 ( − 1)) and in the ray approximation, without allowance for refraction, is described by the integral: where ( , , ) is a function describing the distribution of the refractive index in a medium. The transformation of the rotation of the coordinate system for different probing angles is taken into account under the integral sign. In this approximation, the delay does not depend on z, but the computation of the delay in accordance with Equations (4)-(5) will have an insignificant dependence on z. In fact, the function Δ ( , , , ) is a linearized back shadow projection of the objects under investigation.
Numerically, the operations in Equations (4)-(5) are implemented using the fast Fourier transform algorithm, which allows us to speed up the solution of the inverse problem many times.
Next, we use the analogy with the method of shadow projections used in X-ray tomography [14,15]. According to the method of shadow projections, it is necessary to summarize Δ ( , , , ) for different angles: where ′( , , ) is the cumulative shadow projection. For the simulated scene, a function, ′( , 0, ), was calculated, which is shown in Figure 6b. In this image, one can already distinguish test objects. It can be seen that an object with a large refractive index has a more intense image (a darker color). To compare our method and the ray approximation shadow projection technique, we have calculated the same image with ray approximation ( Figure  6a). Visually, our approach gives a better quality. We performed data processing with additive noise of 30% and 100% (Figure 4b,c). As a result of data processing, images of given inhomogeneities in noise conditions were obtained (Figure 7a,b). One can see that with an additive noise of 30%, the recovery result is little different from the restored image without noise. With an additive noise of 100%, a significant distortion of the image occurs; however, it is still possible to distinguish the presence of simulated objects. Moreover, an object with a large refractive index is visualized more clearly.
We normalize the function ′( , 0, ), believing that ′( , 0, ) = 0 at a refractive index = 1 and maximum values are taken at = 2. Figure 6 shows graphs of the cross-sections of functions ′( , 0, ) for the first and second objects. In the graph (Figure 8), the values of the normalized function ′( , 0, ) approximately correspond to a given distribution of the refractive index.
Furthermore, according to the method of shadow projections for X-ray radiation, it is necessary to calculate the inverse convolution with Tikhonov regularization from a function ( , , ) to a function of the form 1 √ + ⁄ in the following way: is the flat spatial spectrum of the reaction to a pointwise heterogeneity in the summation of the shadow projections. The regularization parameter is selected interactively based on the noise level and artifacts in the restored image; the larger this parameter, the lower the noise level but the more blurry the image is. Value ( , , ) is proportional to the values of the refractive index in the medium. As a result of applying Equation (8), the image of dielectric heterogeneities was calculated ( Figure 9). In this case, the regularization parameter, , is chosen to be equal to 3 ⋅ 10 ⋅ , where is the maximum modulus of the square of the function ( , ).   The refraction coefficient could be calculated by the following formula: ( , 0, ) = ( , , ) The application of Equation (8) allows an increase in the sharpness of the image; however, since the occurring wave effects are not equivalent to the ray approximation as in X-ray tomography, it is necessary to set the regularization parameter to be relatively large, which leads to the blurring of the resulting image. Equations (7)- (8) are an approximate solution of the inverse problem of restoring the distribution of dielectric heterogeneities. In fact, radio waves of backward propagation are not direct rays, as in the case of X-ray radiation. In reality, there are effects of refraction, reflection, diffraction, and multiple scattering. However, in approximating that backward propagation waves are direct rays, Equations (7)-(8) should ensure the reconstruction of images of radio-permeable dielectric heterogeneities. Nevertheless, Equation (4) takes into account the wave nature of the probe radiation.
In addition, numerical modeling was carried out for objects of various shapes. The results of image reconstruction using Equation (7) are shown in Figure 10. The dashed lines in Figure 8 show a simulated objects contour. Figure 10a shows the result of the image reconstruction of two cylinders with refractive indices of 2 and 1.5. Because of the diffraction of waves on one object before falling onto another, the shape of the objects being reconstructed is distorted. As the radius of dielectric cylinders decreases, the diffraction effect no longer substantially distorts the images of objects (Figure 10b). The numerical simulation of a hollow cylinder was performed (Figure 10c), which showed that the circular shape of the object is restored correctly but that the internal volume is restored with distortions, which is explained by the refraction effects at the cylinder boundary. An object of a complex step-shaped form was also considered (Figure 10d). In the area of shading itself, a distortion of the propagation paths of the waves arises, because of which the steps are not clearly drawn.
Algorithm 1 shows the steps for reconstructing the real part of the refractive index.

Algorithm 1
1. The spatial amplitude-phase distribution of the field is restored for all sounding frequencies and for all sounding angles. The reconstructed field is multiplied by the conjugate field of direct wave propagation in a homogeneous medium.
2. The time domain is transitioned to by integrating the reconstructed field in frequency. 3. The time delay of the equivalent impulse relative to the field in a homogeneous medium is determined at each point in space for each angle of sounding by searching for the maximum signal in time.
4. The resulting time delays are considered as shadow projections of objects with different refractive indices. The higher the refractive index, the greater the time delay. In this sense, the sought value is associated with shadow projections in the same way as in the classical tomographic method for X-ray radiation. The shadow projections for various angles are summed up, and a three-dimensional image of the inhomogeneities is obtained. However, this image is blurred by the hardware function of the system. 5. Deconvolution is applied with regularization to take into account the hardware function of the system and a three-dimensional image of inhomogeneities in time delays is restored. 6. A three-dimensional image for converting time delays into the values of the refractive index is normalized.

Conclusions
The authors developed a method for the broadband multi-angle transmission radiowave tomography of dielectric objects based on the backward wave propagation model. A comparison is made between the time delay of the signal calculated from the measurement results in the entire sounding region and the propagation time in a homogeneous medium.
Studies were conducted on the modeling of this method at various objects. On the basis of the results of these numerical simulations, it can be concluded that the method allows the obtaining of a resolution in the order of a wavelength in the theoretical limit. The combination of time delays obtained at different probing angles makes it possible to reconstruct the refractive index distribution in the medium. However, eight projections are sufficient to distinguish between the dielectric permittivities of several objects.
Note that the proposed method does not take into account the refraction of waves in dielectric objects and multiple reflections in the probed region. However, as shown by the results of the numerical experiments, this method is highly resistant to additive noise. Even when adding noise with a level of 100%, it is still possible to visualize objects in the study area.
In the future, the results of numerical simulations will be checked by the results of field experiments. However, even now, we can say that the proposed method and the studies carried out are of great interest for radio tomography.
The novelty of the proposed method consists of taking into account the wave propagation of the probe radiation, and the time delays acquired by the wave when passing through dielectric inhomogeneities throughout the probed area are used as shadow projections.
The proposed method is high-speed because of its non-iterative approach. Good performance is realized through the use of the fast Fourier transform algorithm. The image of heterogeneities is calculated without multiple iterations. This method is not iterative, but we assume that in the future, this method can be used to obtain the first approximation for iterative methods to reduce computation time, but this is a separate study and requires additional research.