Underwater Gas Leak Quantification by Convolutional Neural Network using Images

1. Introduction

Drilling and exploration of wells in the subsea environment make use of pipelines to transport oil and gas. Due to corrosion, structural failure, or third-party damage, the subsea pipelines are subject to material leaks [1,2]. The risks increase especially if the pipeline is in decommissioning state and/or after a long time of service. The consequences can be significant to the ecosystems, which also incur financial losses for all stakeholders [3]. Therefore, monitoring tools are an essential part of a safety culture during the process operation. Several methods have been developed for detecting leakages in different domains [4]: acoustic [5,6,7], chemical [8], pressure [9], and others. Some of these sensors can be attached to remotely operated vehicles (ROV), introduced internally via pigging/crawling robots, or installed on seabeds.

An ROV equipped with an onboard camera can provide real-time filming for monitoring purposes during the inspection of a leak suspicion or routine surveying throughout the production. Since the operation of an ROV can be costly, a fixed camera could be set permanently in strategic positions for continuous surveillance [10]. Nonetheless, in both cases, the images just provide information on whether there is a leak. One could be interested in extracting more information regarding the nature of the leakage, e.g., flow rate or type of fluid leaked. This type of information will influence whether there is a need to intervene in the well. A successful intervention demands trained personnel and financial resources. Hence, for a better assessment, it is of great interest to increase the ROV’s capabilities by developing a quantitative tool to support the intervention decision-making process.

In this sense, Computer Vision is a field dedicated to the interpretation of images by a computer that leads to actions or inferences [11]. Image processing techniques can aid the Oil & Gas industry in solving problems in various applications on refineries and upstream plants [12,13,14,15]. With respect to bubbly flow, Chen et al. [16] developed an algorithm for the processing of bubbles originated from high-speed image experiments and Generative Adversarial Networks (GANs); they coupled background subtraction, binarization (using Otsu’s method), boundary extraction, and arc segmentation.

Regarding subsea leaks, Wang and Socolofsky [17] developed a stereoscopic image system coupled to an ROV to study natural seeps’ bubbles and droplets in deep-sea waters of the Gulf of Mexico. In this case, the laboratory apparatus consisted of air bubbles released in a 1 m tap water tank at 20°C and atmospheric pressure. According to the authors, the air is considered a good representation of methane fluid dynamics for image testing under environmental conditions. The image analysis sequence that the authors performed depended upon the bubble flow type. In short, they employed the Sobel gradient filter, canny edge detection, thresholding, and watershed algorithm. Important leak features were calculated after all bubbles were identified, such as the bubble size distribution and flow rate. The volumetric gas flow rate was calculated based on the average total volume of bubbles in a frame and the average bubble displacement time in the vertical direction. After an initial application of this method in natural seeps in the Gulf of Mexico, Wang et al. [18] and Razaz et al. [19] expanded this method to other observations in the same ocean basin. In another study regarding natural seeps, but in the South China Sea, Di et al. [20] calculated the bubble diameter distribution and release rate employing a semi-automated algorithm that relied on manually determined thresholds. She et al. [21] proposed the gas flow quantification on the marine environment with the aid of median filters, canny edge detection and ellipse fitting. However, as Liu et al. [22] stated, this traditional image processing approach cannot deal with the variety and complexity of different scenarios.

Therefore, despite the widespread applications of the traditional image processing approaches, the Computer Vision field can also rely on new tools arising from Artificial Intelligence (AI), which are learning-driven by data accumulation and have a superior capacity to meet the stringent demands (accuracy, application transfer, efficiency) [22]. Convolutional Neural Networks (CNNs), as a subarea of AI, gained notorious attention in the scientific community for this purpose [23,24,25], since they can efficiently handle data organized in tensors, which includes images [26,27,28,29], time-domain information [30,31,32], radars [33], pressure signal [34], acoustic information collected from sonars [7], and flow information [35,36], together with Computational Fluid Dynamics (CFD) [37,38,39,40,41]. Furthermore, due to the scarcity of data or model complexity, transfer learning may be employed for time-saving and model enhancement purposes, i.e., training a CNN architecture on pre-trained weights. Transfer learning was performed on other chemical engineering applications, for instance, in fault detection [42,43], quantum chemistry [44] and in thermodynamics [45]. Another common question concerns the architecture hyperparameters, which are defined a priori of the training. Neural networks are adaptable to handle various types of data. However, tuning and optimizing different hyperparameters, which is essential for achieving optimal performance, can be time-intensive. The goal of hyperparameter optimization methods is to minimize the generalization error predicted by a machine learning model across different potential configurations. Sadoune et al. [46] combined Bayesian optimization with a tree-structured Parzen estimator to forecast biogas production using real data from wastewater treatment plants.

Hence, CNNs can be applied to study fluid flow, which includes flow regime identification, detection of two-phase flow, and indication of leaks. Concerning identifying bubbles in a water column, Caldas et al. [39] employed CNNs for detecting bubbles of simulated data. The algorithm consisted of a classification problem, yielding 99% of accuracy. Poletaev et al. [35] applied CNNs to identify patterns in a two-phase bubbly jet. They generated 30,000 synthetic images comprising a broad range of overlapping, blurriness, and volume gas fractions. Their methodology is described as follows: (1) classify the probability of bubble presence in the image after comparing 60 CNNs’ different configurations; (2) calculate bubble center employing an algorithm; (3) separate bubbles from the noise background using a denoising autoencoder (a multilayer perceptron network), which facilitates the determination of geometric parameters. The NNs were validated against experimental data. Ilonen et al. [36] compared several bubble-detecting methods, among them CNNs (whose accuracy was superior to the others) in an oxygen delignification process for a pulp production line. The authors also described a bubble size estimator that does not require prior bubble detection. Khaydarov et al. [47] suggested that the bubble size distribution in a stirred bioreactor can be a target variable for enabling smart equipment because it can indicate mixing quality and gas holdup, thus saving operational time. For this purpose, they used a You Only Look Once (YOLO) network to detect air bubbles. The data could be used to complete a CFD model. Another work that used a YOLO was made by Xiang et al. [48], who applied it as a bubble detector in artificially generated bubbly flows. Rutkowski et al. [49] also carried out image processing by the YOLO network and Faster R-CNN in microfluidic droplets.

Some studies were limited to detecting the plume as a whole, not individual bubbles. Shi et al. [26] also applied Faster R-CNN to leak detection in ethane cracker plants. Similarly, for subsea gas leak monitoring purposes, YOLO and Faster R-CNN networks were employed by Zhu et al. [50] and Hu et al. [51].

Nevertheless, all cited applications focused solely on object detection, which in general: (1) demands intensive human labor on manually labeling the dataset during the training phase; (2) it is a task only to predict if some object is present in the image, creating a bounding box close to the object; (3) consequently, it does not dictate the precise object boundaries, requiring further work on region delimitation and contour definition for precise bubble size calculation. Therefore, a more simple and promising alternative is to use semantic segmentation instead of object detection. Semantic segmentation is a classification process on the pixel level, assigning each one a label. Thus, it predicts if some pixel belongs to the category, e.g., to an object or the background. It can be a strategy for multi-phase flow processes in chemical engineering that require knowledge of particle diameter distribution. One of the CNNs developed for this task is the U-Net created by Ronneberger et al. [52] for cell identification in biomedical applications, a field where Computer Vision is well-established. Other networks commonly used for semantic segmentation are DeepLab series, SegNet, and PSPNet [22]. Schäfer et al. [53] employed U-Net to segment oil droplets dispersed in a water column. The post-processing phase enabled a straightforward calculation of the particle size distribution. Furthermore, Caldas et al. [40] employed U-Net for measuring bubble diameter using synthetic images generated by Computational Fluid Dynamics in different conditions of simulated leaks. Cerqueira et al. [54] generated a Particle Tracking Velocimetry (PTV) algorithm using U-Net for identifying and highlighting drop contours of oil in water on a centrifugal pump impeller. Bergau et al. [55] proposed to quantify small methane leaks employing a 3D CNN model, but using artificial laser spectroscopy, not real images.

In this work, we propose an innovative methodology that contributes to the literature as follows:

• to demonstrate an easy way to automatically annotate a huge dataset (77,676 images) for CNN training purposes;
• to show how to perform hyperparameter optimization on U-Net in Big Data situations;
• to apply the U-Net network for gas flow rate quantification and bubble diameter calculation in a water column employing a laboratory model of subsea leaks. In this step, a comparison between transfer learning and hyperparameter optimization is carried out. In the former, the network architecture and the weights from the step based on the synthetic images serve as the starting point. In the latter, a hyperparameter optimization on the U-Net is developed to find the best architecture without any prior knowledge. In this step, the training is performed in a Big Data context. The results are validated against the experimental values obtained in the laboratory.

2. Theory

This section aims to briefly overview the techniques used in this work and already cited in Introduction. First, we will shed light on more traditional image techniques in Image processing techniques. Then, Machine Learning, and more specifically, convolutional neural networks, will be discussed on ML.

2.1. Image Processing Techniques

2.1.1. Thresholding

Thresholding is a segmentation technique in which the classification is based on the histogram of a property of an image, which is usually its grayscale. It aims to find a value that can divide the image into a group with less than this value and another one with greater or equal to this value. It can be viewed as a statistical decision process in which the probability of belonging to a specific class is searched. Hence, the histogram of the image property is treated as an approximation of its probability density function [11].

Suppose an image whose histogram of grayscale intensity $f(x,y)$ presents a bimodal distribution. An example is an image with bright objects on a black background. Therefore, extracting objects in this image is straightforward by finding the value T that divides the pixels into two modes, and consequently, the background and foreground are found. An intuitive approach is to determine a threshold T whose output $g(x,y)$ classifies into [11]:

$$g(x,y)=\left\{\begin{array}{cc}1\hfill & \mathrm{if}f(x,y)\ge T\hfill \\ 0\hfill & \mathrm{if}f(x,y)<T\hfill \end{array}\right.$$

(1)

The output $g(x,y)$ is called a binarized image whose background is denoted as zero, and the foreground is marked as one. If T is a single value valid for the entire image, thresholding is denoted as global. If T depends on the position $(x,y)$, selected independently for a group of pixels, they are called local (or dynamic/adaptive) thresholding [11,56]. This work will discuss only the main global thresholding method: Otsu’s. This method aims to maximize the interclass variance present in the histogram. In this sense, the threshold T should give the best outcome possible between classes regarding grayscale intensity values [11]. The method is described below as introduced by Otsu [57].

Supposing that a threshold $T\left(k\right)=k$ ($0<k<L-1$) is selected so that the image is divided into two classes, $C_1$ and $C_2$; $C_1$ consists of pixels with levels $[0,k]$, and $C_2$ consists of pixels with level $[k,L-1]$. Using the threshold, the probability of each class occurrence are ${\beta }_1$ and ${\beta }_2$, respectively; the mean values of each class are calculated by ${\mu }_1$ and ${\mu }_2$, respectively; and class variances are ${\sigma }_1^2$ and ${\sigma }_2^2$, respectively. In this sense, ${\sigma }_T^2$ is the global variance, and ${\mu }_T^2$ is the mean value for the whole image histogram. Adopting a discriminant criterion $\eta ={\sigma }_B^2/{\sigma }_T^2$ that evaluates the goodness of the threshold at the level k, where the subscripts $BC$ and $WC$ represent the between-class and within-class variances, respectively:

$${\sigma }_{WC}^2={\beta }_1{\sigma }_1^2+{\beta }_2{\sigma }_2^2$$

(2a)

$${\sigma }_{BC}^2={\beta }_1{({\mu }_1-{\mu }_T)}^2+{\beta }_2{({\mu }_2-{\mu }_T)}^2$$

(2b)

$${\sigma }_T^2={\sigma }_{WC}^2+{\sigma }_{BC}^2$$

(2c)

Therefore, the method aims to maximize the between-class variance ${\sigma }_{BC}^2$, which is a measure of dissociation between classes. In this sense, ${\sigma }_1^2$ and ${\sigma }_2^2$ should be the farthest from each other. Since, ${\sigma }_T$ is a constant (${\sigma }_T>0$), $\eta$ is also maximized. Thus, the algorithm should search for the optimum threshold $k^{\ast }$, such that [11,57]:

$$\eta ={\displaystyle \frac{{\left[{\mu }_T\beta \left(k\right)-\mu \left(k\right)\right]}^2}{\beta \left(k\right)[1-\beta (k\left)\right]}}$$

(3a)

$${\sigma }_{BC}^2\left(k^{\ast }\right)=\underset{0<k<L-1}{max}{\sigma }_{BC}^2\left(k\right)$$

(3b)

2.1.2. Sobel Gradient Filter

The main use of the Sobel gradient filter is to find edges on an image by computing their derivatives. The first-order derivatives indicate an edge because maximum or minimum points of intensity represent the sudden change of pixels. The gradient vector of an image reads:

$$\nabla f(x,y)=\left[\begin{array}{c}{G}_x\\ G_y\end{array}\right]=\left[\begin{array}{c}{\displaystyle \frac{\partial f}{\partial x}}(x,y)\\ {\displaystyle \frac{\partial f}{\partial y}}(x,y)\end{array}\right]$$

(4a)

The magnitude of the gradient vector of an image is called the gradient image:

$$M(x,y)=\sqrt{{G_x}^2+{G_y}^2}$$

(4b)

Since images have a discrete representation, it is not possible to proceed as a continuous differentiation. Therefore, the Sobel operators try to approximate the computation of the derivatives by convolving the Sobel kernel with the image, represented as a two-dimensional array A (* is the convolution operator):

$$G_x=\left[\begin{array}{ccc}+1& 0& -1\\ +2& 0& -2\\ +1& 0& -1\end{array}\right]\ast A$$

(4c)

$$G_y=\left[\begin{array}{ccc}+1& +2& +1\\ 0& 0& 0\\ -1& -2& -1\end{array}\right]*A$$

(4d)

Other types of gradient operators are Roberts and Prewitt [11,56].

2.1.3. Fast Fourier Transform

The Fourier transform $\mathfrak{F}$ (FT) for an image $f(x,y)$ with height H and width W is given by Equation (5). The Fourier transform computes the magnitude and phase shift response at each frequency $\xi$ and repeats this procedure all over the domain [58].

$$\begin{array}{c}\hfill \mathfrak{F}\left\{f(x,y)\right\}=F(\omega ,\upsilon )={\displaystyle \frac{1}{HW}}\sum _{y=0}^{H-1}\sum _{x=0}^{W-1}f(x,y)e^{-j2\pi ({\displaystyle \frac{\omega x}{W}}+{\displaystyle \frac{\upsilon y}{H}})},\\ \hfill x=0,1,2...W-1,y=0,1,2...H-1.\end{array}$$

(5)

$\omega$ and $\upsilon$ are the frequency indices in the horizontal and vertical directions, respectively. The component $F(0,0)$ represents the average intensity of the images, located at the center of the frequency domain representation. Higher frequencies are located away from the center, near the edges. The discrete Fourier transform can be represented in polar coordinates as:

$$\mathfrak{F}\left\{f(x,y)\right\}=\left|F(\omega ,\upsilon )\right|e^{j\varphi (\omega ,\upsilon )}$$

(6a)

where $\left|F\right(\omega ,\upsilon \left)\right|$ is the Fourier spectrum and $\varphi (\omega ,\upsilon )$ the phase frequency defined by:

$$\left|F(\omega ,\upsilon )\right|=\sqrt{R{e}^2(\omega ,\upsilon )+I{m}^2(\omega ,\upsilon )}$$

(6b)

$$\varphi (\omega ,\upsilon )=ta{n}^{-1}\left[{\displaystyle \frac{Im(\omega ,\upsilon )}{Re(\omega ,\upsilon )}}\right]$$

(6c)

The use of Fourier Transforms for practical purposes is often infeasible due to computational costs of $O{\left(HW\right)}^2$. To address this issue, the Fast Fourier Transform (FFT) algorithm was discovered, reducing the operations to $HW{log}_2\left(HW\right)$, which led to the development of image processing as known today. One fundamental property of FT that made it possible to implement is the separability of the 2-D FTs. For more details of the algorithm, please consult Gonzalez and Woods [11].

2.2. Machine Learning Techniques

Artificial Intelligence (AI) is the capacity of a computer system to display cognition. AI research has focused on achieving a certain goal by a computer with the best possible outcome, similar to what is done in optimization theory, which minimizes a cost function. The AI subfield dedicated to designing systems that learn based on experience - without being explicitly programmed - is called Machine Learning (ML) [59]. The strategy lies in discovering underlying relationships based on data acquisition.

Artificial neural networks (ANNs) are part of the ML field. A recurrent definition of ANNs found in the literature is that the networks mimic the functioning of the biological brain, in which the connections between neurons - axons and dendrites - are arranged in a massively parallel way and stimulated by synaptic weights, which constitute the learning driving force [60]. In a multilayered ANN, the information flows through connections between inputs and outputs, which are the nodes (neurons) of the ANN organized in layers.

2.2.1. Convolutional Neural Networks

A type of neural network is the Convolutional Neural Network (CNN). CNN became one of the most popular networks for deep learning. Its popularity is due to its ability to process the data in the form of tensors (multidimensional arrays), such as signals and language (1-D), images and audio (2-D), and videos (3-D) [25].

Figure 1 shows an overview of the architecture. The classical convolutional neural networks [24] comprise convolutional layers, pooling layers, and a fully-connected layer. It is a type of feed-forward neural network. The convolution layer extracts important characteristics from input into a feature map by using a kernel (or filter). In a convolutional layer l, the input I from layer $(l-1)$ is processed by several kernels K of size $k\times k$. A single convolution operation at the point $(u+i,v+j)$ involves multiplication of $I_{u+i,v+j}$ by $K(i,j)$ as stated in Equation (7) [61,62]:

$$O_{u,v}^{\left(l\right)}=\theta \left(\sum _{i=0}^{k-1}\sum _{j=0}^{k-1}{I}_{u+i,v+j}^{\left(l-1\right)}·K\left(i,j\right)+b^{\left(l\right)}\right)0<u\le H,0<v\le W$$

(7)

where H and W are the input height and width, respectively; b is the bias; $\theta$ is an activation function, and $O_{u,v}$ is the output. The feature map size is reduced by using a smaller kernel than matrix [61], leading to sparse connectivity, which helps save considerable computational resources. In turn, traditional multilayer perceptrons go through a full matrix multiplication operation between patterns and weights, which could slow down the training process. The number of kernels and their shape are hyperparameters, and each kernel element is a neuron activated by a function defined by the user.

The subsequent layer of a CNN evaluates nearby feature maps by using a summary statistic function of them (average or maximum); this is referred to as pooling. As an effect, dimensionality data is reduced, lessening possible overfitting and memory requirements. It does a type of downsampling operation, reducing the size and aggregating the information forwarded in the subsequent layers. This property can reduce computational costs. They are commonly followed by a dropout layer that randomly deactivates neurons to prevent overfitting. Several convolutional and pooling layers can be arranged together to extract features.

These first CNNs were developed in a sense to be invariant to location, i.e., if a motif appears in a particular part of an image, it could appear anywhere [25]. Since the focus was on pattern recognition, this approach can save computational costs and speed up the classification process. However, in image segmentation, more specifically in semantic segmentation, in which the classification is done pixel-wise, the concern is spatial localization. In order to capture both context and spatial information, Ronneberger et al. [52] developed the U-Net architecture, composed of two parts: a contracting and an expanding path. The contracting path is responsible for the image downsampling, i.e., the feature map sizes are successfully reduced in the subsequent blocks at the cost of expanding the number of filters. The contracting path is organized in blocks of convolutional and pooling layers, which captures the context (for instance, “do we have a bubble present here?”). After that, it is then followed by blocks of deconvolutional layers that perform the inverse convolutional operation and more convolutional layers to maintain the context obtained. These blocks and the skipped connections form the expanding path that reconstructs the feature map to its original size and regains space information. Thus, the expanding part carries out an upsampling process. The skipped connections are copy layers transmitting data from the downsampling section to a twin block in the upsampling one [52,63,64]. At the end of the network, there is another convolutional layer whose output gives a probability tensor of the same image dimension. The value given by a sigmoid activation function (or a softmax if multilabel segmentation) predicts the likelihood of belonging to some class. The U-Net structure is shown in Figure 2.

Batch normalization is a generalized technique applied in the Machine Learning community to improve training. According to the original paper [65], batch normalization aims to stabilize the input changes introduced by updating the preceding layers’ parameters, also known as internal covariant shift (ICS). In each layer input, the ICS is minimized by setting the batch distributions’ two moments (mean and standard deviation) to one and zero, respectively. In turn, Santurkar et al. [66] argue that batch normalization smooths the optimization landscape by making the gradients more predictable and well-behaved.

Another common strategy employed in ML is dropout. It aims to reduce overfitting by randomly deactivating neurons and their connections, temporarily removing them from the network and making the neural network thinner. According to Srivastava et al. [67], this feature enables a neural network to be more robust by preventing co-dependence on nearby connections and to be more exposed to a diversified type of input.

2.2.2. Hyperparameter Optimization

A topic of discussion in ANNs is network hyperparameter tuning. Hyperparameters are the parameters defined a priori from training (in contrast with the weights, which are optimized during training), and they directly affect the performance and the accuracy of the training results. Therefore, the hyperparameter optimization (HPO) problem is composed of several hyperparameters, such as neural network’s architecture, learning rate, batch size, and momentum [68,69]. Yet, it can be a complex problem to solve due to the costly computation (for large/complex datasets), complex configuration space (comprised of continuous, categorical, and conditional hyperparameters), lack of smoothness, and a loss gradient function [70].

Several methods were developed to tackle this issue. One of them is grid search. Grid search is the most basic HPO approach, in which the complete combination of hyperparameters is trained, and the best one is elected based on criteria. However, it can be computationally unfeasible due to the number of combinations to be evaluated. It suffers from the curse of dimensionality. i.e., the number of objective functions to be computed increases exponentially with the search space, which leads to an unfeasible amount of data available for a reliable result. A second approach is to use random search, i.e., choosing randomly hyperparameters based on their importance according to a budget B [69,70]. According to Bergstra and Bengio [71], random search has a low effective dimensionality because the HPO objective function is more sensitive to changes in some dimensions than others.

Bayesian optimization (BO) is an HPO framework that evaluates a black-box function (model-based) to find an optimum set of hyperparameter [69,70]. The optimization goal can be described as searching for the global maxima at the sampling point of an unknown objective function $\psi$:

$${\mathbf{x}}^{*}=\underset{\mathbf{x}\in \chi }{arg\; max}\psi \left(\mathbf{x}\right)$$

(8)

where $\chi$ represents the search space, and ${\mathbf{x}}^{*}$ the vector space of hyperparameter inputs part of $\chi$. Bayesian optimization is based on Bayes’ theorem, which states that the posterior probability $P\left(M\right|E)$ of a model M given an event E is proportional to the prior probability $P\left(M\right)$ and the likelihood $P\left(E\right|M)$ of E given M:

$$P\left(M\right|E)\propto P(E\left|M\right)P\left(M\right)$$

(9)

which is also applicable to BO: it employs some prior information of the function $\psi$ and sample points to evaluate the posterior of $\psi$. This posterior information guides where $\psi \left(x\right)$ is maximized, driving the optimization next step based on a criterion, an acquisition function. Thus, Bayesian optimization consists of two components: a surrogate model fitting all observations to the unknown function, and an acquisition function that decides the region to evaluate next, determining the utility of each candidate, taking into account the exploration and exploitation trade-off [69,70,72].

The most common surrogate model employs Gaussian kernels as prior distribution because of their simplicity and their capabilities to deal with uncertainty. For the acquisition function, the expected improvement function is usually employed [46,69,70,72]. Alternatively to BO, Tree-Parzen Estimator (TPE) models $P\left(M\right)$ and $P\left(E\right|M)$ directly [73,74].

A major disadvantage of BO is that it does not stop unpromising trials from occurring, which, in turn, can be costly from a computational point of view. In this regard, Hyperband is an extension of the successive halving method (and a variation of random search), which is focused on allocating a budget B among the N hyperparameter configurations. First, it evaluates all configurations with a $B/N$ budget for each, then discards the worst half and allocates more resources to promising candidates until one single configuration is left with the maximum budget. Since Bayesian optimization considers prior experiments, and hyperband considers early stopping to save computational time, recently Falkner et al. [75] integrated both Bayesian optimization and bandit-based methods into a single one: BOHB.

2.2.3. Transfer Learning

Transfer learning (TL) is a valuable technique for transferring knowledge across domains, particularly useful when dealing with small datasets, complex hyperparameter tuning, and time-consuming training processes. TL involves using a pre-trained model, one that has already been trained on a large dataset, and adapting it to a new, but related, task. This approach is made by transferring the knowledge captured in the model’s weights and biases to the new dataset, which helps in generalizing predictions effectively [43].

In chemical engineering, TL can be beneficial for tasks such as scaling up processes or modeling similar types of equipment. For example, a model trained on data from one type of reactor can be adapted to predict performance in a different but similar reactor type, leveraging previously learned patterns to improve accuracy and reduce training time [76].

3. Materials & Methods

3.1. Experimental Procedure

A laboratory or a reduced model is a prototype that reproduces the system or phenomena on a minor scale in order to test a hypothesis or a scientific proposal. Nascimento et al. [77] designed an experimental apparatus to generate a reduced model of the subsea leakage, as shown in Figure 3. In this case, the gas employed is air for safety purposes. The air is blown through an orifice in a tap water glass reservoir. The blower consists of an air compressor connected to a voltage regulator and stabilizer. The output voltage sign is responsible for varying the air flow rate. The air flows from the compressor through a flow rate stabilizer and a silicone tube, and then it is injected into the water tank through an orifice. The orifice of injection has different hole diameters $d_i$: 0.5 mm, 1.0 mm, 2.0 mm, and 5.0 mm. The temperature and pressure are on room conditions (20 °C and 1 atm). In the Supplementary Material, Table S1 details the different conditions: the orifice diameters $d_i$, the different types of orifice from where the bubble is released, the flow rate Q set, and the two-phase fluid system (tap water/air and saltwater/air). The apparatus is described below:

• Voltage regulator 115 V manufactured by Varivolt: maximum 11 electric current 11 A, maximum power consumption 1.5 kVA, input voltage 115 V, electric tension variation 0-130 V;
• Voltage stabilizer manufactured by SMS, electric power 300 VA, input current 2.5 A, input and output voltage 115 V;
• Air compressor (positive displacement by diaphragm type) manufactured by Big Air, model A32, flow rate 3.5 L/min;
• flow rate stabilizer manufactured by Cole Palmer, model GV-32908-73, flow rate range 0.5 to 50 liters per minute;
• Glass tank with 96 liters in volume. Dimensions: 40x60x40 cm.

The wall tank was coated with ethylene-vinyl acetate. According to Nascimento et al. [77], the intent was to reduce reflectivity by using a black matte color. A light source was provided inside the tank next to the orifice. The images were captured using a camera from a Xiaomi MI 9T cellphone at a rate of 240 frames per second (fps).

The volumetric flow rate Q is obtained by counting the bubbles collected in a beaker during a specific time $t_c$ recorded by a stopwatch, Equation (10).

$$Q={\displaystyle \frac{Vo{l}_{beaker}}{\mathrm{\Delta }{t}_c}}$$

(10)

According to [77], a possible approach to estimate the bubble diameter $d_b$ is based on the number of bubbles $n_b$ collected within the time interval $t_c$ in the beaker. In light of this approach, the volume of bubbles is calculated considering the bubbles ideally occupy all the volume beaker $Vo{l}_{beaker}$ without empty spaces and overlaps, Equation (11a). The estimation also treats bubbles as perfect spheres in all dimensions, Equation (11b).

$$Vo{l}_b={\displaystyle \frac{Vo{l}_{beaker}}{n_b}}$$

(11a)

$$d_b=\sqrt[3]{{\displaystyle \frac{6Vo{l}_b}{\pi }}}$$

(11b)

3.2. Ground Truth Dataset Generation

Figure 4 shows examples of frames (for an extract of a frame for all 44 video IDs, please consult Fig. S1 in the Supplementary Material), each one representing a video ID. In the laboratory, the experimental images of the reduced model cannot be uniform in illumination, as seen in Figure 4, with either bright or dark spots present in the background. Furthermore, the light intensity does not precisely mark the liquid and gas interface. The background can sometimes be brighter than the bubble, as seen in video #8 in Figure 4a. On the contrary, from videos #10 to #33, the bubble region is clearly brighter than the background (see #14 and #19 in Figure 4b, Figure 4c). In a third scenario, the background and foreground share similar pixel intensities (#34 to #44, see #37 in Figure 4d). Also, another challenge is the presence of other objects that make the segmentation less straightforward, as the back bright “line” seen in #1 - #9 movies. Since the lighting differs not only across the images but especially between the different leak states (represented by a video ID), this suggests that a general method that comprises all these conditions must be developed to carry out segmentation properly. Convolutional neural networks are good candidates for this task.

High-quality data has to be fed for proper training of ANNs, the so-called supervised learning process. This data used as standard is otherwise known as “target” or “ground truth data”. When using synthetic images, the pixel intensity is well defined as described by Caldas et al. [40]. Thus, a simple Otsu’s thresholding method (please refer to thresholding) satisfies as an image segmentation technique to generate the masks for training. Nonetheless, different procedures may be employed to tackle mask generation due to the challenges presented in real images. At first, we proceed with the Global Otsu methodology as described in Caldas et al. [40]. It did not yield good results (an improper value of threshold suggested) as shown in Figure 5 due to the poor light conditions in the area analyzed. For instance, Figure 5a shows that the method could capture part of the information of the bubble contour. However, the penumbra on the left side of the image made it impossible to establish a value that could properly delimit the background and the bubble. It should be noticed that the thresholding was inverted since yellow points represent pixels with the value number “one”. On the other hand, Figure 5b shows better thresholding results, but still, the contours were not adequately defined. Figure 5c depicts a poor segmentation due to a key light difference in the background, which led the algorithm to compute these differences in pixels as the threshold. Figure 5d presents another case of failed segmentation; in this case, because the left side of the image was brighter than the right one; therefore, the algorithm determined this bright side as belonging to the bubble. Evidently, an alternative procedure must be undertaken to overcome the issues mentioned above. A certain degree of proper segmentation when using Otsu’s thresholding suggests that a miscellaneous of techniques may be the solution, which needs to be developed.

Combining the image processing methods described in Image processing techniques was the solution to overcome the pixel intensity heterogeneity. As the first step, the Sobel filter was used to delimit the contours, in which clear contours were found.

Based on these contours, Otsu’s thresholding was used to partition the regions. Then, the morphological operation of closing is used for filling holes and improving contours. Figure 6a and Figure 7 depict the methodology employed and the results in this case, respectively.

Nevertheless, the “Sobel/Otsu/Closing” algorithm returned poor results for videos #17, #25, #33. The reason for that can be found after analyzing the flow conditions: leaks that flow under intermittent regimes, i.e., bubbles that do not flow continuously, led to poor identification of contours by the Sobel filter since some frames captured were empty, without any bubbles. Thus, another high-pass filter should be employed. In the videos above, the Fast Fourier Transform was used to highlight the boundaries. When creating the “FFT or Sobel/Otsu/Closing” algorithm, the consideration used for choosing either the FFT or Sobel filter on each frame was based on whether the image was blank. A decision criterion was developed based on a value T, manually tested for each of the videos. Figure 6b illustrates the flowchart and Figure 8 or Sobel/Otsu/Closing shows the final results from the “FFT or Sobel/Otsu/Closing” method. It should be noticed that the method was capable not only of segmenting the bubble and the background but also of returning an empty image when no flow was present, as depicted in Figure 8c.

Table S2 in the Supplementary Material summarizes the methods employed for each video. It was decided not to use videos #40 - #44 because the camera was not static, which would compromise the measurement of the equivalent bubble diameter. Other videos which could not produce good true masks were discarded (01, 04, 07, 38, 39). For converting the area computed as a bubble to a physical scale, the height of the image was measured by a ruler and compared to the height of the frame in pixels, which resulted in a conversion factor $fc$.

3.3. Training CNNs with Real Images

In this part of the work, two methodologies are compared employing U-Net. The steps followed for each one are:

• Importing the U-Net architecture and weights from a synthetic image model performed by Caldas et al. [40]. No hyperparameter optimization is carried out. This step is referred to as “transfer learning”.
• Building new U-Net models based on hyperparameter optimization of the architecture and training, i.e., number of filters on a block, activation functions, whether to apply batch normalization, dropout value, optimization method, learning rate, and weights decay. The weights are trained from scratch, i.e., no prior information is given to the network. This step is referred to as “hyperparameter optimization”.

In method 1, the network structure is maintained. A total of 7,759,521 parameters are trained. The first block of the contracting path has 32 filters of size 3x3, and the number of filters is doubled on each block, reaching 512 filters in the most profound block. Then, the number of filters is halved in the expanding path. On each block of the contracting path, a convolutional layer employing ReLU as an activation function and a dropout of 20% is succeeded by a max pooling layer and by another convolutional layer with ReLU. The contracting path is composed of blocks of transposed convolutional layers and convolutional layers. The last convolutional layer has a sigmoid as an activation function. The batch size is 32.

In method 2, the hyperparameter optimization employed a combination of Tree-Parzen Estimator (TPE) to provide sampling and a Hyperband algorithm as the watchdog for early stopping. The Optuna framework was the library used [78]. The search space is described in Table 1. A total of 432 configurations are possible. However, the idea is not to perform a grid search but to obtain a fast response with fewer experiments. One-fifth of the dataset was employed for time-saving purposes during the hyperparameter optimization. Then, the best configuration found was trained again on the whole dataset.

In both cases, the dataset split was 60%, 24%, and 16% for training, validation, and testing.

As a metric to train the neural network, we employed the Dice-Sørensen coefficient, Equation (12a), a measure that computes the spatial overlap between two segmentation sets. Representing $y_n$ as the ground-truth sample n, $\widehat{y_n}$ the corresponding predicted sample, and $i,j$ indexes relative to points in an image, the following expression can be written:

$$DSC(y_n,\widehat{y_n})=2{\displaystyle \frac{|y_n\cap \widehat{y_n}|}{|y_n|+|\widehat{y_n}|}}={\displaystyle \frac{2{\sum }_{j=0}^{256}{\sum }_{i=0}^{256}{\widehat{y}}_{ij,n}·y_{ij,n}}{{\sum }_{j=0}^{256}{\sum }_{i=0}^{256}{\widehat{y}}_{ij,n}+{\sum }_{j=0}^{256}{\sum }_{i=0}^{256}{y}_{ij,n}}}$$

(12a)

For a dataset with N samples, the global Dice-Sørensen coefficient is written as in Equation (12b).

$$GDSC(y,\widehat{y})={\displaystyle \frac{2|y\cap \widehat{y}|}{\left|y\right|+|\widehat{y}|}}={\displaystyle \frac{2{\sum }_{n=1}^N{\sum }_{j=0}^{256}{\sum }_{i=0}^{256}{\widehat{y}}_{ij,n}·y_{ij,n}}{{\sum }_{n=1}^N{\sum }_{j=0}^{256}{\sum }_{i=0}^{256}{\widehat{y}}_{ij,n}+{\sum }_{n=1}^N{\sum }_{j=0}^{256}{\sum }_{i=0}^{256}{y}_{ij,n}}}$$

(12b)

It should be highlighted that the loss function considered for training was $1-DSC(y_n,\widehat{y_n})$, in which each sample and mask are computed individually. Then, the average for all samples and masks is calculated [79]. This strategy, rather than using the global metric as loss — Equation (12b) — was used to force the network to learn each image individually, not the whole batch. Since the background composes most part of the data, it has a significant impact on the loss calculation. Thus, missing one bubble in the whole batch may have a negligible effect on the global Dice, whereas in an individual image, missing a bubble may be noticed by the individual Dice [79]. It is a workaround to avoid the well-known problem of class imbalance [80].

It should be stressed that the U-Net is incapable of quantifying by itself. The neural network only supplies a predicted mask with defined contours. The bubble diameter ${\widehat{d}}_b$ is computed in a post-processing step. Yet, it is straightforward since regions that belong to the bubble have a pixel with code “one”, and the ones to the background have a pixel number “zero”. For this purpose, the Regionprops function was employed in the Scikit-image library [81]. This widespread function is also present in other imaging libraries, like OpenCV and MATLAB. This function calculates the area based on the pixels, as explained before, and the equivalent diameter is based on the area:

$${\widehat{d}}_b=\sqrt[2]{{\displaystyle \frac{4A_b}{\pi }}}$$

(13)

For measuring the flow rate, first, the major axis a and minor axis b of the bubble are measured in a determined horizontal position. Since the area of the bubble can vary during upwelling, it was chosen to use two positions ($y=120$ px and y=220 px from the top). The flow meter only calculates the geometric properties if a bubble passes through. Since the bubble is an ellipsoid, the volume is calculated as follows [20]:

$$\widehat{r}=\sqrt[3]{{\displaystyle \frac{a^{2}b}{8}}}$$

(14a)

$${\widehat{V}}_{eq}={\displaystyle \frac{4}{3}}\pi {\widehat{r}}^3$$

(14b)

When the bubble passes through the flow meter, the average of the volumes is considered to avoid multiple counting.

3.4. Computational Resources

All code was written in Python 3. For generating the masks described in ground-truth, the library Scikit-image was employed [81]. The neural network models were built through the high-level deep learning API Keras [82] in the Tensorflow v2 platform [83]. For the hyperparameter optimization, as mentioned above, we used the Optuna framework [78], including their TPE and Hyperband algorithm implementations. All experiments were monitored with Keras callbacks and, when applicable, Optuna callbacks.

We carried out the training and hyperparameter optimization in a Dell XPS 8940 built with an Intel^® Core^™ i7-10700 CPU, 16 GB RAM, 8 cores; and coupled with a GPU NVIDIA^® GeForce^® RTX^™ 3060 with 12 GB of memory.

4. Results & Discussion

In this section, the numerical flow rates and bubble diameter values obtained by transfer learning and hyperparameter optimization are presented. For the validation of the model, the equivalent bubble diameter ${\widehat{d}}_b$ is compared with the experimental ones, $d_b$, done by Nascimento et al. [77] as shown in Equation (11b). In this section, the analysis is separated into groups, which are organized according to the type of flow, orifice, and fluid. The first has a continuous flow, capillary tube as orifice type, and tap water as fluid (Table 2, where $d_i$, Q, $v_o$ are the orifice diameter, flow rate, and release velocity, respectively). The second one has the same conditions but has an intermittent flow behavior. Group 3 has other types of orifice but has a continuous flow, and the system is composed of tap water (Table 3). The last one, group 4, has salt water, continuous flow, and a capillary tube as an orifice type (Table 4). Only in groups 1 and 3 were the experimental bubble diameters $d_b$ available.

Furthermore, the bubble diameter values predicted by the literature in the same conditions are presented, i.e., with the same leak flow rate and orifice diameter. It should be highlighted that these equations are only valid for the first bubble released in a water column and not for the average of bubbles in a water column. They were described here just for qualitative purposes, not for validation. Next, the probability distribution function is computed based on a specific statistic value for the bubble diameter. Finally, some computational considerations are also discussed in this work. It should be noted that, during the discussion of the results, the frames were shuffled during the split of training, validation, and test datasets (the training of the network does not consider the temporal relationship). This step is necessary to guarantee a good fit for the CNN training. Therefore, the results are not presented in the same order as the original videos. Furthermore, the dataset split was done to avoid neural network overfitting, as explained in meth:cnn. Since the training dataset has 60% of the original dataset, the results will be discussed using this set, except where noted.

4.1. Computational Performance

Figure 9 the history of the transfer learning training. Regarding computational performances, a quick convergence was found in terms of loss. The algorithm stopped in less than 25 epochs by employing the early-stop approach. Furthermore, the training loss ($1-DSC(y_n,\widehat{y_n})$) achieved a value of less than 0.10 in the second epoch, which indicates that the training is efficient due to the large data available. The average of the Dice-Sørensen coefficient computed for the image and targets individually reported were 0.91 and 0.68 for training and validation, respectively, while the global Dice-Sørensen coefficient was 0.97 for training and validation. The difference in global and local metrics can be explained due to the major presence of the background (class imbalance in semantic segmentation), as explained in meth:cnn.

Figure 10 illustrates the hyperparameter optimization and the values achieved in each trial. The pruning algorithm allows a quick choice of promising candidates. We ran 64 trials that lasted for 110.2 h. Out of 64 trials, 41 were pruned, and 23 were complete. The early stopping works two-fold: by using the pruning algorithm itself and by using the callback criteria of stopping if no improvement $1\times {10}^{-4}$ is achieved during 3 epochs. Therefore, only 7 trials (#2, #9, #21, #47, #54, #56, and #58) lasted more than 10 epochs of training. Trial #58 lasted 20 epochs, which also traduced in more computational time (6.62 h) and had the best validation loss.

Figure 11 illustrates the parallel coordinate plot. The training returned, in general, elevated values of the target metric ($>0.7$). It can also be seen that the combination of Adamax as an optimizer, dropout of 0.3, and “False” on dropout on upsampling is a bad choice. All three candidates that contained this combination figured in the six worst ones in terms of metric performance. To investigate further, Table 5 shows the 10 trials that demonstrated the best performance for the global Sørensen-Dice coefficient. It can be inferred that using batch normalization is vital for high performance. A good balance between computational cost and precision is needed to find an optimal model. In this sense, the number of filters substantially influences the time necessary for computing as the number of weights to be trained increases. Therefore, the best architectures were selected to be trained through the complete dataset.

To better understand the influence of each parameter, Figure 12 shows the parameter importance (in terms of the fraction of variance explained) using a functional Analysis of Variance (ANOVA) evaluation [84]. As expected, batch normalization exerts the most substantial influence among all parameters, followed by the activation function, though the former has a much greater effect. On the other hand, the number of filters has a minor influence, as evidenced by the contour plot of the hyperparameter combinations shown in Figure 13. In this figure, the number of filter charts is, in most cases, symmetric, whereas when batch normalization is used, better results are seen. The optimizer method selected still has some effect: studies that employed “Adam” showed poorer outcomes depending on the combination. For dropout, architectures using 0.1 and 0.2 obtained, in general, superior performance. With respect to the learning rate, both 0.01 and 0.001 values shared similar results. Using dropout on upsampling seems to be more promising than not using it. Regarding the activation function, the “selu” function may be the best choice. Therefore, based on the discussion above, an apparent successful combination matches the architecture of trials #58, #56, #55, which are the first, sixth, and eighth best hyperparameter performances.

In deep learning training, the learning rate directly influences the generalization of the model, which trades off the overfitting and underfitting [68]. As a rule, a minimal learning rate can lead to overfitting and slow performance. In our model, the learning rate range was stable and had little influence over the metric. For instance, for the configuration found in trial # 9, changing the learning rate from 0.001 to 0.01 impacted the metric from 0.8780 to 0.8758. However, the computational cost is significant; the first lasted for 3.96 h, while the latter lasted for 0.84 h. This result is according to the findings in the work of Smith [68], who studied the influence of several hyperparameters in benchmarking datasets. The author suggests that an optimal learning rate can speed up the training without becoming unstable. The results of the complete training from hyperparameter optimization are discussed in the following sections. In general, they were very similar to the results of the transfer learning step.

4.2. Group 1

The numerical values of bubble diameter ${\widehat{d}}_b$ obtained by hyperparameter optimization are shown in Figure 14. A good agreement is found between the median values and the experimental ones, found by counting the bubbles in a beaker [77]. Using the maximum bubble could provide some overestimation for the bubble diameter. The exception is on the largest orifice diameter ($d_i=5.0$ mm) and shown in Figure 14c. However, it should be noted that the experimental values for this case also lie within the median diameter error bars for 3 out of 7.

Furthermore, after plotting the values predicted by the correlations from the literature [85,86,87], the median values were found either above them (but close to the experimental values, Figure 14a and Figure 14b), or their error bars lie in the region of the literature curve (Figure 14c). Moreover, the values are consistent for each CNN dataset — training, validation, and test — which shows the algorithm’s robustness. Therefore, the methodology presented herein is considered valid for computing the bubble diameter.

The bubble diameter curves from the literature [85,86,87] reported here are for the detachment from the nozzle. The bubble diameter in a column can change due to the force dynamics. An increase in bubble diameter along the column can occur due to the decrease in static pressure (considering no breakup/coalescence) [88]. It is a hypothesis for overestimating the values from the expected in the literature.

Figure 15 shows examples of prediction masks and the corresponding bubbles printed. By analyzing this figure, it can be inferred that the CNN trained by transfer learning could reproduce the ground truth data (or the true mask) well. In addition, in Figure 15a, the U-Net did not consider other objects (see the left bottom corner) as bubbles, which is a positive performance. The algorithm is capable of processing multiple bubbles successfully, even for high-speed leakages, like in Figure 15a. However, it should be noted that dense bubbly flows with several overlapping bubbles were not considered in this work.

The bubble probability density function (PDF) was estimated using kernel distribution [89]:

$${\widehat{f}}_h\left(x\right)={\displaystyle \frac{1}{nh}}\sum _{i=1}^{n}K\left({\displaystyle \frac{x-X_i}{h}}\right)$$

(15)

where n is the number of samples, h is a smoothing parameter known as bandwidth, and K is the kernel used — for Gaussian: $K(x;h)\propto exp(-{\displaystyle \frac{x^2}{2h^2}})$.

The PDFs for the training dataset for median values of the bubble equivalent diameter are shown in Figure 16. The median values were chosen here as they are the most representative statistic based on the preceding discussion. After analyzing Figure 16a and Figure 16b, they present a high probability density (all with $>0.10$) with a well-defined bell shape, which suggests a low level of randomness. However, for the largest orifice, as shown in Figure 16c, the bubble diameter median values have a low density, which indicates an elevated level of randomness, as already pointed out by the greater error bars.

Figure 17 describes the median bubble equivalent diameter values on each frame for the different orifice diameters. There is a clear trend towards more variability in larger orifice diameters, as noted in Lu et al. [90], who studied the influence of orifice size diameter on an external loop airlift reactor. The authors also observed a decrease in the probability density by increasing the orifice diameter. The smaller peaks in a PDF, as seen in Figure 16, are equivalent to a less uniform bubble diameter.

Transfer learning yielded results similar to those of hyperparameter optimization. For instance, the bubble diameters predicted only diverge in the second decimal number, and the figures generated are almost identical to Figures 14-17. Please refer to Figs. S2 to S5, which shows the comparison of bubble diameter values with the experimental ones, the average of the median values of bubble diameter, the diameter of the bubble predicted, and the PDFs for transfer learning.

4.3. Group 2

Figure 18 depicts the bubble diameter median values for the group with intermittent flow for hyperparameter optimization. The CNN could correctly assign images without bubbles, which resulted in several images with diameters’ values equal to zero. Furthermore, it correctly captured the frames when the bubble was present. It should be highlighted that these videos contained few bubble releases. Thus, the algorithm recognized a sequence of no-flow frames, which resulted in zero bubble diameter. Since the dataset was split randomly, the Figure 18 was also with random values. The average values of the bubble diameter, not considering the frames without bubbles, were 3.93 mm for video #33 ($d_i=1.0$ mm), 3.36 mm for video #17 ($d_i=2.0$ mm), and 6.73 mm for video #25 ($d_i=5.0$ mm). It should be noted that only one bubble is released in this intermittent flow regime. Therefore, a comparison with the literature on bubble detachment is valid. The model of Jamialahmadi et al. [86] predicted 3.5 mm, 4.3 mm, and 5.5 mm, respectively, while equations from Gaddis and Vogelpohl [87] indicated 3.5 mm, 4.5 mm, and 6.0 mm, respectively. The underestimation in the video from $d_i=2.0$ mm can be attributed to the slower release from the nozzle, in which the U-Net calculates the diameter even before the bubble detaches from the nozzle.

Transfer learning presented some differences from hyperparameter optimization. The average for the videos, not considering the frames without bubbles, were 3.97 mm for video #33 ($d_i=1.0$ mm), 3.37 mm for video #17 ($d_i=2.0$ mm), and 1.64 mm for video #25 ($d_i=5.0$ mm). Again, the underestimation in video #17 is because the bubble is formed and stays on the nozzle for some time before total release. In video #25, small objects ($<2.0$ mm) contaminated the water column and were confounded with bubbles by transfer learning. Results from video #25 for hyperparameter optimization are much closer to the ones expected by the literature because the small objects were correctly not considered as bubbles. The median values for group 2 using transfer learning are shown in Fig. S6.

4.4. Group 3

The bubble diameter results from group 3 after trial #58 are displayed in Table 6. The results showed that $d_i=1.0$ mm was again the one that found the closest match to the experimental values, with less than 2% of deviation from the average of median predicted bubble diameter values. For $d_i=5.0$, the algorithm underestimated ${\widehat{d}}_b$, but all values lie within the standard deviation, which was larger in both cases. It should be highlighted that the values expected from the literature were divergent for $d_i=1.0$ and $d_i=5.0$. The correlations of Akita and Yoshida [85], Jamialahmadi et al. [86], Gaddis and Vogelpohl [87] expected 3.5-3.7 mm for video 02, 3.9-4.6 for video 03, 6.3-7.0 for video 08, and 7.0-9.5 for video 09. Hence, it could be inferred that the measurement of bubble diameter for larger orifice diameters is an object of divergence. Especially, Akita and Yoshida [85] does not agree with other authors for higher flow rates.

Figure 19 illustrates the predicted masks for group 3. This group is different from the other group in terms of the proximity scale with the bubble. The camera is much closer to the bubble flows. The U-Net was capable of reproducing the true mask well. For instance, in Figure 19a, the line behind the bubble was correctly considered as background in the image. This consistency of accurately measuring bubbles from different scales and perspectives demonstrates the potential and flexibility of the algorithm.

Comparing both approaches in Table 6, it can be seen that the results are similar for $d_i=1.0$ and videos 03 and 09, but closer to the experimental values for videos 02 and 08 for the hyperparameter optimization. The standard deviations stay in the same numerical range for all cases analyzed. Therefore, the CNN could differentiate again very well the foreground (bubble) from the background even when the true mask incorrectly assigned some pixels to the foreground (please refer to Fig. S7, which shows predictions for transfer learning in Group 3 for other conditions than in fig:HOpredictedexamplesGroup3). The reason for that is the large amount of information the CNN was exposed to, which led to a precise recognition of objects.

4.5. Group 4

Regarding bubble equivalent diameter, Figure 20 shows the probability distribution function and median values for Group 4 for hyperparameter optimization (for transfer learning, please refer to Fig. S8). Table 7 compares the predicted values of bubble diameter with values from literature [86]. The results display very low deviations from the equivalent diameter as indicated by the high density in the PDF, as seen in Figure 20b. Of all the groups studied, Group 4 had the lowest deviation. Figure 21 shows the different output masks predicted by transfer learning (for a similar figure by HPO, please refer to Fig. S9). The contours were again well delimited. However, a small object on the left in Figure 21b is noticed. Despite being capable of differentiating some objects as background, circular objects’ pixels could be misclassified as bubbles, which is a limitation of semantic segmentation. An object detection algorithm must be coupled with the U-Net to avoid this.

4.6. Flow Rate Prediction

Considering the fact that, in a frame, some bubbles can be in the time step of entering that particular frame, the algorithm can underestimate the bubble diameter for that specific bubble. Therefore, the maximum and median bubble diameter values of a frame were chosen as the representative statistic of this frame. Then, the average of the bubble diameter of each statistic value was computed in the different conditions of flow rate Q and orifice diameter $d_i$.

The predicted flow rate by training CNN by transfer learning against the experimental ones is illustrated in Figure 22a. The values obtained by the algorithm correspond closely to the observed ones. Out of 35 points, only 8 are above the 20% deviation line.

As mentioned, trial # 58 was selected for training through the whole dataset for the hyperparameter optimization. The results of the predicted flow rate against the expected one are shown in fig:HO-flowrate-predicted. A good accuracy is found between them, with a coefficient of determination of 0.941.

5. Conclusion

In this work, we presented a methodology for measuring bubble diameter for subsea leaks. Besides that, we added several novelties to the literature, such as the automatic generation of ground truth data and hyperparameter optimization for U-Net. The results obtained are in line with the data obtained in the literature both by transfer learning and hyperparameter optimization. The algorithm uses the binarized image not only to calculate the area but also to generate outlines of where the bubble is present, along with a screen print of the diameter. The good prediction capacity of gas flow rates in the early stages can significantly enhance process safety. According to Figueredo et al. [3], most subsea accidents involve the release of gas fluid, which leads to loss of lives, explosions, fire, and environmental impact. In the case of Brazil, the authors pointed out the lack of integrity assessment by the operating companies after analyzing National Agency of Petroleum, Natural Gas, and Biofuels (ANP) reports. The present study helps to tackle this issue by a direct evaluation of gas leakages from images in the initial stage (minor leaks), which is a more proactive approach than remediation after propagation. Thus, the methodology is applicable to leaks in deep waters, and the application of the algorithm in the Oil & Gas Industry is immediate.

Furthermore, considering the transfer learning and the hyperparameter optimization, both approaches yielded similar results with a slightly superior performance for the latter. Hyperparameter optimization could be more computationally intensive, but it seeks the best possible outcome in terms of topology design. On the other hand, transfer learning is a potential route when a numerical simulation is available.