1. Introduction
Industrial tomography (also known as process tomography) allows the following process to happen inside the industrial tank or pipes without direct contact with the process itself. The tomographic problem is part of a family of problems known as inverse problems. The main goal of tomography is to find a solution to the tomography problem, which gives the image of the examined object given the measurements collected at the edge. Thus, tomography enables monitoring and a better understanding of industrial processes and facilitates real-time process control. Thanks to developed, non-invasive techniques, it is possible to study the flow of fluids and the density of substances in each area in two or three dimensions. In tomographic systems, spatial resolution is the ability of the instrument to distinguish various objects of different shapes and sizes. In tomography, a better resolution is achieved with the right combination of several sensors, reconstruction method, and model geometry [
1].
Typical applications of process tomography in the industry include monitoring various concentration profiles in mixing and separation vessels in single-phase and multiphase processes. Process tomography allows to control and understand processes in real-time. In addition, tomography provides current feedback on processes, their effectiveness, and response progress. As a result, it can reduce the amount of waste and improve the company’s overall energy efficiency. The model of analysis and control of industrial processes is shown in
Figure 1.
In the crystallisation processes, changing the reaction conditions can lead to a significant increase in production capacity but requires a significant change in many aspects of process design and consideration of unforeseen difficulties or benefits. An important aspect is to characterise the individual stages of the process. Process tomography sensors can provide quantitative and qualitative measurements of the dynamics and volume of underlying processes. Data can be used as basic information about the process. Tomographic sensors can also help in determining the mixing efficiency and other performance characteristics. Tomography can be used to continuously monitor process reliability by providing information as to when reaction conditions change. The installation of a two-phase flow identification system can be extremely useful in controlling the presence of air bubbles in liquids that are intermediates, in which air is undesirable.
In many cases, the presence of air bubbles in production processes can cause irreversible losses. Additionally, the air in liquid intermediates with higher viscosities in the pharmaceutical, chemical, and food industries can be very unfavourable in some cases. Real-time control to detect the presence of air bubbles is therefore necessary. Industrial tomography is used in processes such as crystallisation and fermentation, which are key chemical reactions in the pharmaceutical, food, and chemical sectors [
2]. Tomographic sensors can explore the entire reactor area and provide information on reaction kinetics and concentrations. It can help indicate the endpoints of the reaction more accurately and provide information on how to ensure the proper chemical conditions for the proper conduct of the process. Tomographic measurements can also be used in various process ranges to confirm scale characteristics.
Currently, a standard for industrial production processes is to monitor the process using sensors that only provide local measurements, which can be insufficient for precise control. Therefore, there is a need for measurements with a high spatial solution. Such processes are designed in block models with the energy and mass exchange elements.
Many parameter codes result in a loss of spatial prediction capability and the status of the current process conditions. Additionally, the accuracy of the reconstruction decreases due to natural physical complexity in phenomena such as fluid dynamics, crystallisation, or fermentation processes. Nevertheless, liquid concentration distribution, phases, and chemicals can be studied. Furthermore, the data obtained can monitor process reactions and improve quality, efficiency, and flow rate [
3].
The automation of industrial processes increases the level of control over their implementation. It guarantees stable, automated, and flexible work. Furthermore, with the technology shown in the paper, one can monitor and obtain measurements at any given moment. Thus, it allows automation and better control of the quality of production processes.
The most popular monitored values for such processes are temperature, pressure, or fluid flow speed. However, the most valuable information obtained is what is happening inside the reaction tanks in chemical industrial processes. This information is needed to control technological processes more efficiently. The only way to achieve this information is with the use of the tomography method. Thus, there is room for the ultrasound tomography method, where a limited number of inexpensive ultrasound transducers are used. The biggest challenge in ultrasound tomography is the reconstruction with the use of the obtained time-dependent data.
Obtaining proper image reconstruction requires that ultrasonic tomography must meet two primary conditions. First, a sufficiently significant difference between the acoustic impedance of different substances has to be obtained to measure an incomplete reflection (a complete reflection will block the possibility to obtain imaging from deeper layers). Second, the wave speed distribution in the measurement area under study should have linear changes since it would be impossible to obtain the depth of these reflections from wave time data with the nonlinear changes.
In industrial applications, tomography-based systems are used mainly for quality control [
4,
5]. Advanced industrial process control analysis can be applied to internal fluid distribution, multiphase flows, fermentation, or crystallisation processes using tomographic sensors. However, industrial tomography applications typically present challenges obtaining reconstruction images, given the process boundary’s observational data [
6].
Ultrasound is a prevalent method used to diagnose and treat patients, mainly in cancer treatment and antenatal diagnosis [
7,
8,
9,
10,
11,
12]. The main problem with today’s ultrasound devices (incredibly portable ones) is the poor quality of the reconstructed image [
13,
14]. There are three main aspects of diagnostic ultrasound physics: spatial resolution, temporal resolution, and contrast resolution. Excessive damping is associated with loss of amplitude for both low and high-frequency transducers used in reflective ultrasonography. Temporal resolution is constrained by the depth of penetration and size of frames, and contrast resolution is constrained by image memory. The paper [
15] shows the physical aspects of potential improvement of reconstructed image quality, but we focus on developing mathematical image reconstruction methods even with cheap transducers.
The existing methods for reflective ultrasound tomography analyse the reflection on the path typical to the transducer [
16]. In contrast, other methods use a convolution-back projection algorithm for reconstruction [
17], or iterative method involving minimisation of the objective function [
18] or use different discretisation model [
19]. The only limit here is the number of transducers used for the ultrasound measurement, and thus it is limited to the amount of information extracted from the system.
In literature, various numerical methods can solve tomography problems for the 2N configuration of transducers. In addition, various information from the signal can be used to solve this problem [
5,
20,
21]. Unfortunately, none of these solutions can give an accurate solution. Most algorithms depend on the transmitted signal from which the amplitude and time of flight (ToF) are used. The problem is solved with the given data algorithms involving the inversion of a non-square matrix, which can only give an approximate solution to the problem [
22,
23].
The machine learning approach is widely used in computer tomography, but most applications are in transmission tomography, such as X-ray or Magnetic Resonance Imaging (MRI), where various Deep Learning types and procedures are used [
24]. Predefined dictionaries of image samples are used in well-known algorithms, such as Global Dictionary Learned in the Statistical Iterative Reconstruction (GDSIR) or Adaptive Dictionary Learned in the Statistical Iterative Reconstruction (ADSIR) and many others [
25], where the image is reconstructed based on low-dose sinogram with different types of Convolutional Neural Networks (CNN) equipped with various types of decoders and encoders. These methods are called image domain reconstruction [
26], where deep learning techniques are only used to improve reconstruction quality and remove artefacts from images. Much less common are data domain methods in which deep imaging is applied to sinograms, and the reconstructions come directly from raw data. There is also a dual-domain method (or hybrid domain), combining image and data domain methods. However, all the above mentioned methods use analytical formulas to recover the image from projection, called filtered back projection (FBP).
This paper shows deterministic and machine learning methods of image reconstruction based on ultrasonography signals. We present a new high-resolution and fast method for ultrasound reflection tomography that is easy to implement, and it is faster because the matrix inversion methods are not needed to obtain the transformation matrix. Furthermore, simultaneously Extreme Gradient Boosting method was applied to image reconstruction as an alternative method. It allows us to boost the reconstruction quality and use the data to visualise the object’s shape.
The outcome of the presented work is a complete system based on the cloud and developed hardware, which allows for constant industrial process monitoring.
The system’s heart is a software platform that integrates elements responsible for the management and communication with the smart devices. The system consists of:
The system is based on cloud architecture and can handle high traffic. Additionally, due to the nature of cloud architecture, heavily loaded elements can be dynamically and automatically scaled to meet the requirements. Part of the system is installed on-premise and communicate and cooperate with part of the platform located in the cloud. The platform is responsible for proper device management and data processing.
The monitoring process using tomography consists of the data collection, transmission, validation, and processing using specialised algorithms.
This paper is organised as follows. The first part of Material and Methods contains a detailed description of the hardware used in the experiments conducted in the laboratory. This section also presents the experimental setup used in the study. The second part of Material and Methods describes the image’s reconstruction idea and starts from the raw signal transformation to the description of two image reconstruction methods based on the second wave-packet reflection time matrix. Both deterministic and machine learning methods are described extensively. Finally, the results section presents reconstructions performed on the same objects by both methods. Any inconveniences and problems associated with the reconstruction of the images are described in the discussion. In conclusion, a summary of the achievements presented in the paper has been presented. Potential spaces for the development of this technology are also indicated.
2. Materials and Methods
2.1. Hardware
The circular tank with the inclusions was measured using the ultrasound tomograph constructed by Research and Development Centre Netrix S.A., Lublin, Poland (3). Measurements were done with the use of the constructed device that can generate and then measure ultrasound signals with the use of replaceable transducers. The idea of the constructed tomograph and the experimental setup is shown in
Figure 2,
Figure 3 and
Figure 4.
The reflective ultrasound tomography is designed in a modular manner. It creates the possibility of easier and faster testing of critical system components. With ready-made modules, we limit the risk of damaging the whole device to its specific modules.
The most important module consists of the main board combined with an analogue signal conditioning board and a display. The second module consists of a specialised high voltage driver, and it is also equipped with a 64-channel multiplexer.
Digital buses of the modules are connected using the RJ45 cat. 6a cables. M12 standard connectors were used, one for every eight channels. The communication with a computer is performed through the USB HS 2.0 standard. The tomograph can save data on external storage devices such as micro-SD cards and USB drives.
The boards are shown in
Figure 5a,b form a multiplexer module equipped with a high voltage symmetrical rectangular signal generator (max about 142Vp-p). The board is controlled by an STM32F103CBT6 micro-controller that is responsible for generating waveforms with high-voltage MOSFET drivers. The software can adjust frequency, the number of pulses generated, and time of shorting the signal output to the ground using the differential CAN 2.0 bus.
A built-in micro-controller does the described analogue multiplexer keys control (CAN bus) or directly by an external microprocessor via an SPI bus. Control buses, i.e., differential “trigger” signal for initiating keying, CAN bus, or SPI bus, is in standard RJ45 cat. 6a (shielded twisted pair). The driver circuit has a 4-channel MD1822 driver chip and double TC8220 bridges. The three-state control allows the transducer oscillations to decay faster after the keying process is complete. The multiplexer module has four basic 16-channel analogue switches for reflection measurements and four additional analogue switches for transmission measurements. Three cascaded hardware counters operating in One Pulse mode were used to generate the appropriate waveforms to control the MOSFET keys.
When the TRIGGER line’s differential input is set to on, the system generates an interrupt in the microcontroller, the ENABLE line is set high, and TIM1 (OUTA, OUTB) is set high triggered. After some time, the TIM1 hardware triggers TIM2, which counts down half the period of the first-timer. Next, the hardware sets low on the OUTA, OUTB lines and triggers TIM3. The third counter counts down when the TX output is shorted to ground and then activates a software interrupt that sets the ENABLE line low. It causes the MD1822 chip to set safe states on the lines controlling the TC8220 high voltage MOSFET keys. The STM32F103 microcontroller, apart from the possibility of cascading counters and generating complementary signals in pulse mode, has the function of setting the dead time—a delay between the opening of the upper and lower key reduces the delay the losses generated during switching. This time it was selected experimentally, and it is not adjustable from the user interface. Another proper function, which found its application in controlling the keys, is the commutation function. Right after starting the counter, it allows programming how its outputs will be set as a commutation event (in the project, the event is the Update Event, i.e., the end of counting the set number of impulses).
Due to the system’s characteristics, we found it is essential to use shorting to the ground method. This way, the generated ultrasound waves have better characteristics, and the reflected wave’s width is smaller than the signal generated without shorting. In the process of design, much emphasis was put on the universality of the main module. The dimensions of the display determined the size of the board as the main module board is mounted on the display. The board has battery backup for the RTC clock, backup registers, and FRAM memory. Communication with a computer can be performed in USB High Speed 2.0 standard. There is also a debugger and serial port connector. CAN, SPI, and RS485 buses are connected via RJ45 cat. 6a connector. Other I/O ports are located on popular pinouts.
The analogue signal conditioning module is an additional expansion board mounted on four spacers. Similarly, a smaller board with a physical layer of ETHERNET communication can be mounted on a smaller terminal strip. For projects requiring a considerable number of I/O ports and where communication via ETHERNET is not required, one can prepare expansion modules clipped onto both terminals. The board can also be used at the stage of prototyping and preparation of other devices. On the board, there is also a small audio amplifier that amplifies the sound from the 8-bit synthesiser built into the FT811 display, which was used to play short characteristic beeps when buttons are touched. The FT811CB display features a 5″ diagonal screen size, 800 × 480 resolution, and a capacitive touch layer with multitouch support. The display’s embedded controller from FTDI features the Embedded Video Engine (EVE2). It allows for efficient display control using only the SPI bus, engaging a small percentage of the microcontroller’s computing power. The analogue module is used as a measuring part of the device. Its role is to obtain a high-resolution signal received by the transducers. The main amplification path was made using the AD8331 measurement amplifier. Gain control at this stage is done by changing the voltage at the GAIN input using a DAC built into the STM32F746 chip. The amplified differential signal is then fed to two separate measurement paths. One of them is a fast LTC2202 10MSps analogue-to-digital converter with a parallel 16-bit output.
In the second measurement path, the signal’s envelope is made by the ADL5511 circuit and is additionally amplified by the AD623 measurement amplifier. The amplification of this path is performed using adjustable digital potentiometer MCP4017T. The circuit converting the signal to an envelope requires an appropriate selection of the value of filtering capacitors FTL1-4. The signal prepared in this way is finally sampled by three analogue-to-digital converters built-in STM32F746 micro-controller, working in triple interleaved mode. It allows the sampling rate of the signal to be increased three times. It results in a sampling rate of up to 7.2MSps (3 × 2.4). However, setting the maximum ADC clock in the micro-controller would require reducing the timing of the entire processor. Thus, the maximum timing for this application (envelope measurement) was 5.1 MSps (3 × 1.7). The timing of both transducers was done using an internal hardware counter. This solution allows for precise adjustment of the sampling rate.
The device is also capable of doing defectoscopy, which is helpful as a standard material analysis tool. The defectoscopy mode is shown in
Figure 6.
2.2. Mathematical Models
The classical tomography methods are based on the idea that the measurement can be translated to the pixels or finite elements via the equation
where
m is the measurement vector,
is the vector corresponding to the pixels or elements, and
is the sensitivity matrix that can translate between both. Unfortunately, the size of
m and
b is not the same in practical applications. Thus, the resulting
A matrix is not square. Finding the inverse of a non-square matrix is a complex computational task. However, some methods can reduce the size of
b [
27,
28]. Although these methods make finding the inverse of
more accessible, they still lack the resolution required to recognise the measured objects’ shape and correct size.
Here we proposed completely different approaches. The deterministic approach uses analytical image reconstruction methods based on the starting arrival time (SAT) matrix. The machine learning approach is also based on the SAT matrix to find reconstruction using the XGBoost algorithm to each finite element of the grid. Machine learning methods used in later method have nothing to do with image filtering or sharpening, but instead tell us whether a mesh point should be marked as a part of the reconstructed object.
Figure 7 presents the idea of how the signal generated by a transducer interacts with the object inside a measurement tank. The path of the first reflection is
I →
R1. The path of the second reflection is
I →
T1 →
R2 →
T2.
Both methods are based on the localisation of wave packets of the raw signal. For each measurement obtained from the tomograph, the SAT matrix is calculated. The first wave-packet is the one that travels in the shortest path (a straight line) from Tx to Rx. The third usually corresponds to the reflection from the inside of the object and should be rejected, as the physical parameters of an object’s interior are generally time-dependent and unknown. The most crucial image reconstruction is the second wave packet, which corresponds to the reflection from the outer boundary of the object, and it is used in the analysis.
On the first stage of the analysis, the raw signal from each probe was transformed into a smoother version of its own by local variance transformation
where
denotes the variance in the neighbourhood of
k-th observation,
denotes the cardinality of the set,
,
mean, respectively,
and
, where
d is the neighbourhood radius of the
k-th observation,
n is the number of analysed time points. Thus, by
we mean the signal value at point
k, and
means the arithmetic mean computed for the neighbourhood of
k-th observations.
The local variance defined in this way reproduces very well the local amplitude variations of the raw signal. The advantage of this solution is that the resulting signal is smoother than the envelope, which translates into an easier finding of the beginning of the reflected wave.
Calibration of this representation can be done by manipulating the size of the neighbourhood radius. To ensure that both the raw signals and their local variances are expressed in the same units, the local standard deviation, the square root of the variance, is used instead of the local variance. As shown in
Figure 8, transforming the raw signal by the local standard deviation did not change the location of the local extrema. The local extrema itself was determined numerically because we do not have an explicit form of the function that determines the local standard deviation. By local maximum
we mean the location where the condition
where
denotes the local standard deviation, and
δ denotes the radius. It should be mentioned that among all the local maxima existing for a given signal, only those that satisfy an additional condition were selected
This condition rejects the local maxima of the signal smaller than a predetermined percentage of the maximum local standard deviation. The solution adopted in this way allows for varying the sensitivity of the proposed tool, depending on the threshold level t. To localise and reconstruct the shape of a phantom, we should find the start of the second wave as a propagation time of the reflected signal. We choose this point by subtracting the radius from the second extremum. In this way, the whole matrix of the reflected signal propagation times is created.
The second stage of the analysis is the reconstruction part, which can be done in two ways. The first one is called deterministic and consists in finding all points (
x,
y) inside the
Ω area for which the following inequality occurs
where
,
are Euclidean distance of (
x,
y) from
,
,
,
are coordinates of transducers,
sound speed of the medium, and
δ is our tolerance.
The method used in the second stage of our analysis was devoted to reconstructing a phantom in the cylinder using a machine learning algorithm. The 2D image of it expressed all experimental setups of phantom localisation and shape. First, we divided our images into 200 times 200-pixel grid, which gives us 28,228 pixels in the cylinder. Then, pixel-by-pixel, we fitted the XGBoost model to predict if this pixel should be marked as a part of the phantom based on the propagation time matrix.
The term “Gradient Boosting” in the Extreme Gradient Boosting (XGBoost) originates from the paper Greedy Function Approximation: A Gradient Boosting Machine [
29,
30]. Unlike random forests, where independent decision trees are a prediction tool in classification and regression tasks, XGBoost models are ensembles of decision trees of a different design.
XGBoost models are widely used by researchers, especially by data scientists, to achieve state-of-the-art on different machine learning challenges. This model performs very well in various prediction and classification tasks, even when the relationships between the dependent variable and the predictors are nonlinear and complex. The mentioned flexibility of the models is because all the boosted tree models are based on an appropriate selection of decision trees and learning them on different data sets. Adding to this, the regularisation of the loss function ensures that these models can describe almost any complex relationship, and there is no overfitting. The biggest disadvantage of this type of model is the lack of explicit form (Black-Box model). However, in the present work, this disadvantage is not particularly important because the analysis of the model’s sensitivity to a change in the value of the SAT matrix is not relevant.
In XGBoost models, each successive tree learns to predict the residual value formed in the previous iterative step. The process of learning the model is based on the minimisation of the objective (loss) function enriched with a part causing the regularisation of the
where
L is the loss function, and
Ω is a regularisation term. For the classification task, the most common choice of
L function is
where
is the observed value of a dependent variable, and
is the predicted value. The regularisation part takes different forms, depending on the constraints we want to apply to the model. These can be regularisation in
L1 and
L2 norms. In most cases, the regularisation part is as follows
where
is the vector of model parameters,
is the prediction of
j-tree, while
γ and
λ are costs of tree depth and model complexity. By setting these parameters appropriately, individual trees are pruned, and the variance of the model is reduced. The learning procedure is usually based on the gradient descent method.
The XGBoost models were built using native libraries implemented into the R language [
31]. All computation and image reconstruction in terms of machine learning was done in the R environment using the following libraries: xgboost [
32], mlr3 [
33], doParallel [
34], tidyverse [
35], plot3D [
36], rio [
37], and splus2R [
38].
4. Discussion
The presented results show that it is possible to obtain high accuracy with few transducers and relatively cheap ultrasound tomography. While the deterministic algorithm can recover most of the details, the machine learning algorithm is an excellent method in finding the location of reconstructed objects.
However, it is not easy to compare the two methods using objective measures of object reconstruction, because the reconstruction proceeded slightly differently in both cases. In the deterministic method, only the edges of the object were reconstructed, whereas in the method using machine learning, both the edges and the interior of the object were reconstructed (this can be observed in the reconstruction plots—see
Figure 10,
Figure 11 and
Figure 12). The differences mentioned above could significantly affect the results of comparisons using measures such as Mean Square Error, Relative Image Error, or Image Correlation Coefficient [
41,
42].
The authors plan to test the developed solution concerning smooth changes in the fluid density. However, this problem is much more complex, and the entire measured waveform will be analysed to perform the reconstruction. It is an exciting topic for further investigation that is out of the scope of this paper. Nevertheless, a review of the results reveals an interesting observation. Regardless of the reconstruction method used in the second stage, its quality depends most on the accuracy of the SAT matrix.
Table 1 describes the results obtained for different methods, including classical ones (Gauss-Newton) and machine learning methods. Of course, not the same quality indicator can be used in all the papers, but we can compare them because they were used in the context of tomography (electrical impedance, capacitance, and ultrasound tomography) because the idea behind them all is to solve the inverse problem.
In
Table 1, the accuracy is the ratio of correctly reconstructed elements on the mesh to the total reconstructed elements on the mesh. Similarly, Relative Error and RMSE are calculated in the context of how close the reconstructed values are to the ground truth. It means that despite different names, all the quality indicators show the “error” concerning ground truth values.
Our method achieved 99.52% accuracy, which can be translated to the value of around 0.48% of Relative Error and 0.0036 RMSE and, as shown, is the best method of the shown methods.
In addition to the results presented in this paper, the authors also performed test reconstructions when the SAT matrix contained only reflection time information for transducers placed close together. Reconstructions based only on information from transducers with small distances between them give equally good results as reconstructions based on the full SAT array. By closely spaced transducers, we mean those whose angle between the transmitting and receiving probes does not exceed 70 degrees. It is most likely because indicating the front of the second wave packet in these cases is much easier than where transducers are placed on opposite sides. However, finding the front of the second wave packet is not easy in that case, and as it is shown that it is not crucial for image reconstruction. A similar solution is used in rotating the ultrasound tomograph, where the transmitting and receiving transducers are located in a rotating transducer [
26].