**3. Results**

#### *3.1. Pool Experiment*

#### 3.1.1. Experimental Design and Equipment

In this study, a scene of leaking gases from a pipeline was simulated in a pool, and the acoustic data was collected by an HT300 PA MBES developed by Harbin Engineering University, China. As shown in Figure 3, the MBES is mainly composed of an underwater subsystem (i.e., the sonar head), an interface unit, and a computer. The underwater subsystem is responsible for the emission, reception, and real-time processing of the underwater acoustic signal. The interface unit is used for real-time information collection of the auxiliary equipment and network data exchange, and the HTCS 3.4 display and control software of this MBES can be operated on a computer. The sonar system emits an acoustic wave with a frequency of 300 kHz, a CW pulse signal with a pulse width of 0.1 ms, and a ping rate of 20 Hz. In addition, its beam width is 2.5◦ × 1.5◦ and max swath width is 126◦.

The distribution of the bubble size and number of bubbles changes continuously during their rising process because they are affected by factors such as water depth, nozzle size, pressure in the pipeline, flow rate, ocean current, and so on. Also, the backscattering strength (or sound attenuation) of the bubbles is closely related to the signal frequency, bubble size, and the number of bubbles [32–35]. For a nozzle, its size decides the sizes of bubbles. For example, the radius of the detached bubble could be expressed as a non-monotonic function of the nozzle diameter [36], and the radius showed an increasing trend during the change of the nozzle diameter (0.1 mm–2 mm). For a single bubble, the different sizes correspond to different resonance frequencies, and the backscattering strength reaches a peak at this frequency [35]. The relationship between the scattering cross-section σ and a bubble of a radius *a* can be expressed as [37]:

$$\sigma(a) = \frac{4\pi a^2}{\left[\left(\frac{f\_0}{f}\right)^2 - 1\right]^2 + (ka)^2} \tag{11}$$

where *f*<sup>0</sup> is resonance frequency of the bubble and *k* is the wave number. This model was applied in [37] to analyze the relationship between backscattering strength and bubble radius at 40 kHz, 180 kHz and 300 kHz. In actual applications, MBESs operate at frequencies from 12 kHz to 400 kHz [3,37], and they have been applied in researching the detection of gas leaks and other targets in water. In general, in shallow water (<200 m), MBESs with a frequency range of 200–400 kHz are usually used for marine surveys, and MBESs with these frequencies are also suitable for the deployment of the unmanned underwater vehicles due to their size. Because the purpose of this paper is to detect underwater pipeline leaks in offshore waters (usually less than 200 m in depth), and comprehensive consideration of factors such as maximum sounding capability, resolution of the multibeam sonar, and frequency response characteristics of the target of interest, this type of MBESs was finally used for experimental research.

During the experiment, the sonar head of the MBES is fixed at 0.5 m below the water surface with a tilt angle of α = 30◦ , as shown in Figure 4. The maximum coverage angle of this MBES is from −63 degrees to +63 degrees and the depth of the pool is shallow; thus, in order to enable the sonar to accurately study the process wherein the bubbles rise to the surface, it is installed vertically with a 30-degree tilt. The leaking gases are generated by an air compressor and released through a long rubber tube, which is connected to the nozzle at the bottom of the pool (5 m depth). Figure 5 is a video capture of the gas leaking scenario during the experiment.

**Figure 3.** A photograph of HT300 PA MBES.

**Figure 4.** Schematic diagram of the pool experiment for simulation of leaking gases.

**Figure 5.** An underwater video screenshot of the bubbles being released.

3.1.2. Optical Flow Calculation of the Water Column Images (WCIs) Containing Leaking Gases

Figure 6a,b are D-T images taken via the conventional beamformer (CBF) and the MVDR beamformer, respectively, and both matrix data are normalized and displayed in decibels. The fixed nozzle diameter is φ = 0.6 mm; the output intensity and the output flow rate of the air compressor are 0.1 MPa and 6.3 L/min, respectively; and the sonar head is stationary. In Figure 6a, fan-shaped bright strings are formed at the MSR, and their brightness is close to the intensity level of gas emissions. In this experiment, the gassy release area is designed outside the MSR and its purpose is to consider whether the detection of the gas emissions is still valid when they occur in an area that have significant interference from sidelobe contamination. Figure 6b shows that the image at the MSR (i.e., the region where the sidelobe effect energy is the largest) was well-suppressed. There also is sidelobe leakage outside of the MSR, but the energy level is usually lower than that of the leaking gases.

Using Figure 6b and its adjacent frame images, the motion field distribution of the obtained image region is an implementation of the Farneback optical flow algorithm in MATLAB with Scheme A. In Figure 7a, the estimated velocity field is displayed in conjunction with the WCI of Figure 6b. The direction of the arrow in the figure is indicated as the direction of velocity, and the length of the arrow is the magnitude of the velocity. Since an estimated velocity value can be obtained for each pixel of the D-T image, to avoid a large number of arrows being too dense to observe, the velocity matrix here is subjected to equal interval sampling processing to display. The image of the small frame on the left side in Figure 7a is the enlarged result of the image in the red frame in the gas release zone. Figure 7b shows the estimation results of the leaking gas motion obtained by processing the same data using Scheme B, and the parameters used by the Farneback algorithm are the same as those in Figure 7a. The estimated velocity is converted to the form of *ux*, *vy* for display. Then, Figure 8 shows a set of WCIs and motion field estimation results after changing the experimental conditions, where the nozzle is replaced by one with a diameter of 1.2 mm, the output intensity of pressure is set to 0.5 MPa, and the output flow rate is 51 L/min.

**Figure 6.** The WCIs containing gas emissions, obtained by the conventional beamformer (CBF) beamformer (**a**) and MVDR beamformer (**b**).

**Figure 7.** The joint display of the WCI and velocity field estimated by frame A (**a**) and frame B (**b**), respectively. The fixed nozzle aperture is φ = 0.6 mm and the output intensity of pressure and the output flow rate was 0.1 MPa and 6.3 L/min, respectively.

**Figure 8.** The joint display of the WCI and velocity field estimated by frame A (**a**) and frame B (**b**), respectively. The fixed nozzle aperture is φ= 1.2 mm and the output intensity of pressure and the output flow rate are 0.5 MPa and 51 L/min, respectively.

It can be seen from Figures 7 and 8 that the velocity field trends estimated by the two processing schemes are roughly similar, but the calculated velocity field matrices have different distributions in the D-T image. The velocity field data calculated by Scheme A and the individual pixels of the D-T image are sometimes in one-to-one correspondence to their spatial position. However, since the velocity field matrix calculated using Scheme B is equally spaced when using the coordinate system with the beam angle and the TWTT as the coordinate axes, the velocity field data near the sonar head is dense, but gradually becomes sparse as the distance increases. In addition, the velocity field of the region containing gas can be visually different from those of other regions. The direction of the velocity in the region containing gas is mostly in the quasi-vertical direction. Since the sonar head is stationary, the ensonified area on the bottom of the pool has not changed. Under this condition, the energy of the individual pixels corresponding to the bottom area in different ping images changes very little; thus, the estimated velocities in this area will be very small, which can be excluded by considering the magnitude of the velocity.

To quantitatively analyze the consistency of the estimated velocities of the two schemes, a square with a side of 5 cm is selected as an observation area to calculate magnitudes of the velocity at different times in this region. The center point coordinate of the observation area is (5.65 m, 3 m). Figure 9 shows the variation of the average of all the estimated velocities in the observation area, along with the sample time. The horizontal axis is the sample time or the image frame number, and the longitudinal axis is the average of the magnitudes of the velocities. For each pixel, the corresponding magnitude of the velocity can be calculated by *V* = *u*2 *<sup>x</sup>* + *v*<sup>2</sup> *<sup>y</sup>*. The curves plotted in Figure 9a,b respectively correspond to the results of the two settings of the nozzle sizes and the output intensities of the pressure conditions, while the data processing is identical. It can be seen from the two figures that the estimated velocity values of the two frames are relatively similar in general, but there are also a small number of different cases. This may be because the number and spatial distribution of velocities estimated by the two frames in the observation area are different. When the difference in velocity values in this area is large, the two statistical averages will show a large deviation. However, this difference does not change the nature of the relatively high-speed movement of the gas emissions from the pipelines, and thus it will not affect the subsequent detection of gas emissions. For leaking gases from undersea pipelines, the velocities of bubbles in water columns were composed of two parts [38]. One is the free rising velocity of the gas bubbles, and the other is the initial velocity when they were ejected at the nozzle. Considering these contributions, the authors of [38] used Doppler shifts of ultrasonic scanning to estimate the leaking gases simulated in the laboratory; the estimated velocity is 3.2–9.7 m/s, which is similar to the results of Figure 9. The echo intensities at different flow rates are further measured and are shown in Figure 10. In the two data acquisitions, the nozzle diameter is always 0.6 mm; the flow rates of the air compressor are 21 L/min and 6.3 L/min, respectively; and the multibeam sonar is set with the same parameters, including the transmit power, fixed gain, and time varying gain. As shown in Figure 10, the large flow rate corresponds to stronger echo intensity as a whole, which is similar to the results found in [37]. This phenomenon implies that there is a certain correlation between the flow rate and the backscattering signal. Compared with results in [39,40], the derived rise velocities in Figure 9 appear on the lower scale. This may be due to the limited output pressure capability of the used air compressor and the absence of currents, which affects the intensity and direction of the motion velocities.

**Figure 9.** The magnitudes of the velocities estimated from Scheme A and Scheme B. (**a**) and (**b**) correspond to the two experimental conditions of Figures 7 and 8, respectively.

**Figure 10.** Comparison of echo intensities at two flow rate outputs.

#### 3.1.3. Detection of Gas Emissions

Due to the WCIs and velocity field distributions in the pool experiment, the thresholds of the image amplitudes and velocity values were set according to Equation (10) to detect gas emissions. The weight coefficient of the amplitude threshold is 1, and the weight coefficient of the velocity intensity threshold is 1. By further eliminating the data below the thresholds in Figures 7 and 8, the results can be seen in Figures 11 and 12, respectively. The velocity field data used in Figure 11a, Figure 12a was estimated using Scheme A, while the velocity field data in Figures 11b and 12b were estimated using Scheme B. As can be seen in Figures 11 and 12, the image amplitude and velocity intensity threshold could be used to exclude images with small velocities but strong backscattering.

**Figure 11.** Results obtained by screening the data in Figure 7 using thresholds of image amplitude and velocity intensity. (**a**) and (**b**) correspond to the processing results of Scheme A and Scheme B, respectively.

**Figure 12.** Results obtained by screening the data in Figure 8 using thresholds of image amplitude and velocity intensity. (**a**) and (**b**) correspond to the processing results of Scheme A and Scheme B, respectively.

#### *3.2. Lake Experiment*

#### 3.2.1. Experimental Scenario

To quantitatively analyze the estimated velocities of the two schemes under the same conditions, the above-mentioned pool experiment is a static measurement where the relative position of the measuring platform and the gas-releasing device is constant. A dynamic experiment is further designed here to measure the simulated gas emissions from a pipeline using an MBES mounted on a survey vessel. The experimental site is in Songhua Lake, Jilin Province, China, and the MBES and air compressor used in the experiment are the same as in the pool experiment. During the experiment, one end of the rubber tube connecting the nozzles is fixed on a triangular holder seated at a flat bottom of water depth of about 13 m, while the other end is extended from the bottom to a barge on the coast and connected to the air compressor. The fixed nozzle aperture is φ= 1 mm and the output intensity of pressure is *Pout* = 0.68 MPa. When the survey ship is traveling, the MBES collects the underwater target echo at a ping rate of 16 Hz.

#### 3.2.2. Detection of Gas Emissions

The WCIs in Figure 13 are obtained when the survey ship travels above the leaking gases, and the velocity field distributions of Figure 13a,b are obtained using Schemes A and B, respectively. From the partial enlargement of the two figures, the corresponding velocity field of the regions containing gas is mainly on the upward trend, and the high magnitudes of velocity are also found in the seabed area. However, as analyzed in Section 2.2.3, the velocity field directions in the seabed area are mostly similar to the seabed terrain tendency. Correspondingly, the velocity direction threshold is set to 40 degrees. For instance, the data in Figure 13a are judged only by the amplitude threshold to obtain the detection result as shown in Figure 14a, and the detection result after additional estimation of the magnitude of the velocity and direction threshold is shown in Figure 14b.

**Figure 13.** The joint display of the WCI and velocity field estimated by frame A (**a**) and frame B (**b**), respectively, from the experimental data of a measurement period in Songhua Lake.

**Figure 14.** The result of the comprehensive detection of the data in Figure 12 using thresholds of pixel intensity, velocity magnitude, and velocity direction. (**a**) is the detection result using only the amplitude threshold, and (**b**) is the detection result of adding the magnitude of the velocity and direction threshold.

It is possible to eliminate pixel data with relatively low intensity, low speed, or high horizontal speeds by comparing Figures 13 and 14 after the detection threshold is set using both the energy and the displacement of the image pixels. Thus, an image with gas emissions can be detected from the complex WCI data.

#### **4. Discussion**

In this paper, the motion field in continuous multi-frame WCIs was estimated using an optical flow method. If only the intensity information of each pixel in the WCI is used to detect the leaking gases, the results could detrimentally include the seabed and other strong interference targets in the water column. Although the images of the seabed and sidelobe interference could simply exclude the portion outside the MSR as was done in [12,15], this strategy may also remove parts of the gassy image, and any strong targets suspended within the MSR region cannot be identified and ignored. As shown in Figure 15, only the amplitude threshold was used to detect the WCI data measured on a short survey line in the Songhua Lake. It can be seen from the figure that it is possible to detect other false targets using only the amplitude threshold, thereby affecting the autonomous and high-accuracy recognition of the gas emissions. For this reason, the estimated motion vector information was further used to eliminate the image pixels of relatively low-speed targets or targets whose backscatter is high but moving horizontally, such as a flat seabed or part of a school of fish. Inevitably, small spots of discrete distributions, as shown in Figure 14, may appear after the processing. These spots are mainly residues after the above threshold decision, and they have a low number of pixels with a discrete distribution. Here, these abnormal pixels can be eliminated by a small sliding window. As shown in Figure 16, the three-dimensional morphology of the leaking gases with a relatively clean background was obtained after multi-threshold detection, which achieves the automatic detection of gas emissions. Considering that sidelobe contamination may adversely affect the amplitude and velocity threshold detection, the MVDR beamforming algorithm was used in the beamforming process to suppress the sidelobes.

In addition, two kinds of schemes were designed for data processing. Similar to [18], in Scheme A, the velocity of the movement is estimated from the D-T image, so that the dimension of the velocity vector matrix is consistent with the D-T image, which is more convenient for visual observation. However, considering that the data processing of Scheme A is mainly a serial structure, and the automatic and real-time detection of the leaking gases on the sonar hardware platform is also a potential demand in the future, this paper further proposes directly estimating the gas motion using T-A images and to judge whether the gas emissions are present or not. In this way, the velocity field estimation can be processed in parallel with the D-T image generation and the seabed terrain detection process from the perspective of the data processing flow, and it is not necessary to wait for the conversion of the image from T-A to D-T, which improves the feasibility of real-time realization in the hardware platform of the MBES.

To verify the effectiveness of the proposed method, experiments were carried out on pools and lakes with relatively shallow water depths. The underwater condition may be more complicated with the increase of water depth. For example, the backscattering strength of the bubbles is also related to water depth (or ambient hydrostatic pressure *P*0). The resonance frequency in Equation (11) can be derived [37] as:

$$f\_0 = \frac{1}{2\pi a} \sqrt{\frac{3\gamma P\_0}{\rho}}\tag{12}$$

where ρ is the density of the water, and γ is the ratio of specific heats. It can be seen from the combination of Equations (11) and (12) that the scattering cross section is also a function of the hydrostatic pressure. This is especially the case in deep water conditions, where the effect of hydrostatic pressure is more obvious. Therefore, when the background noise level is isotropic and the leakage rates and sonar frequencies are the same, the signal-to-noise ratio corresponding to the echo of gases at different depths will also change. At this point, the amplitude threshold of Equation (10) should vary with the depth. Therefore, to fully test the detection performance and adjust its adaptability with changes in depth, in future research we will further carry out experimental verification studies in deep water environments. Of course, the rising gases may be deflected by ocean currents [4], and be affected by water density, salinity, temperature, and static pressure, which are important for seismic oceanography [41]. Therefore, the impact mechanism should be further analyzed and the velocity direction threshold adjusted according to the above conditions in sea experiments.

In a pool experiment, the MBES is stationary. When experimenting on the lake, the MBES is installed on the survey vessel and moves along the survey line. With fluctuation of the vessel, the sonar produces changes in orientation, such as roll and pitch, which will bias the estimated target motion displacement. Therefore, when processing the data on the lake, motion attitude information is compensated. In addition, in the pool experiment, the position of the illuminated bottom is unchanged; thus, the fluctuation of the echo from bottom in WCIs of the adjacent pings is very small. When a survey vessel is sailing, the fluctuation shows more significant changes, which is mainly due to statistical fluctuation of the seabed echo energy. Therefore, in terms of the bottom part of the velocity field in Figures 8 and 13, the estimated velocity in the lake experiment appears to be larger than that in the pool experiment.

This paper aims to detect moving targets through high and low-speed fields, without considering the rigorous verification of the accuracy of the estimated speed. This will be considered in subsequent studies, such as using video-in-situ observations [13] for comparison. We will also consider the combination of the backscattering strength of the bubble group to quantitatively estimate the flow rate of the leaking gases. In addition, it is assumed that the seabed terrain is relatively flat when the velocity direction threshold is used. If the terrain changes steeply, the appropriate values of this threshold may be difficult to obtain. Moreover, fish may be one of the major reasons for misdetections, although the

direction threshold is set. Therefore, in future research, seabed detection technology will be applied to judge the seabed and eliminate seabed image interference. In addition, other features, such as target shapes, will be considered to eliminate interference from fish.

**Figure 15.** The result of the amplitude threshold detection of WCI data measured on a short survey line in the Songhua Lake.

**Figure 16.** The result of the multi-threshold detection of WCI data measured on a short survey line in the Songhua Lake.

#### **5. Conclusions**

To achieve automatic detection of leaking gases from underwater pipelines, a detection method based on the combination of motion and intensity information of WCI pixels was studied in this paper. The motion of the image pixel was estimated using the Farneback optical flow principle, and two data processing schemes were designed according to two types of WCI data. Through the data analysis of two simulation experiments in a pool and a lake, it can be seen that the velocities obtained based on the above two schemes had relatively good consistency. In addition, after the comprehensive threshold

selection using amplitude and velocity information, leaking gases could be detected more efficiently. From the perspective of our designed structure, Scheme B is more suitable for real-time implementation of sonar hardware platforms, while Scheme A is more suitable for execution in display and control software, and in post-processing.

In this paper, it is assumed that the trend of seabed topography is relatively flat in the design of the velocity direction threshold. In future research, bottom-tracking technology will be introduced to judge and eliminate the seabed image with more complex trends of topography, and morphological features will also be introduced to judge and eliminate the other high-velocity motion and strong-backscatter targets in water columns. Moreover, the influence of sound velocity changes for the thresholds will be analyzed.

**Author Contributions:** C.X., T.Z., and W.Z. developed and designed the experiments; C.X., W.Z., M.W., and W.D. performed the experiments; C.X. and W.Z. analyzed the data; C.X., T.Z., and M.W. wrote the paper, and J.L. and P.R.W. modified and polished the paper. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was supported by the National Natural Science Foundation of China (grant numbers U1709203, 41606115 and 41876100), the National Science and Technology Major Project of China (grant number 2016ZX05057005), the National Key R&D Program of China (grant number 2016YFC1402303), and financial assistance from the Postdoctoral Scientific Research Developmental Fund of Heilongjiang (grant number LBH-Q18042). The work of Jianghui Li and Paul R. White was partly supported by the European Union Horizon 2020 research and innovation programme under grant agreement number 654462 (STEMM-CCS).

**Acknowledgments:** The authors would like to thank Yongxiang Hu and Dongyang Li for their tireless work during the experiment on the lake.

**Conflicts of Interest:** The authors declare no conflict of interest.
