Extraction of Submarine Gas Plume Based on Multibeam Water Column Point Cloud Model

Abstract

The gas plume is a direct manifestation of sea cold seep and one of the most significant symbol indicators of the presence of gas hydrate reservoirs. The multibeam water column (MWC) data can be used to extract and identify the gas plume efficiently and accurately. The current research methods mostly start from the perspective of image theory, which cannot identify the three-dimensional (3D) spatial structure features of gas plumes, reducing the efficiency and accuracy of detection. Therefore, this paper proposes a method for identifying and extracting the gas plume based on an MWC point cloud model, which calculates the spatially resolved homing of MWC data and constructs a 3D point cloud model of MWC containing acoustic reflection intensity information. It first performs noise suppression of the 3D point cloud of the MWC based on the symmetric subtraction and Otsu algorithm by leveraging the noise distribution of the MWC and the reflection intensity characteristics of the gas plume. Then, it extracts the point cloud clusters containing the gas plume based on Density-Based Spatial Clustering of Applications with Noise (DBSCAN) according to the density difference between the gas plume point cloud and the background MWC point cloud and next identifies the point cloud clusters by feature matching based on fast point feature histograms (FPFHs). Finally, it extracts the gas plume point cloud set in the MWC. As evidenced by the MWC data collected from gas hydrate enrichment zones in the Gulf of Mexico, the location of gas plume extracted by this method is highly consistent with that of gas leakage points measured during the cruise. Using this method, we obtained the point cloud data set of gas plume for the first time and accurately characterized the 3D spatial morphology of the subsea gas plume, providing technical support for gas hydrate exploration, subsea gas seepage area delineation, and subsea seepage gas flux estimation.

1. Introduction

Seabed gas seepage is a natural phenomenon widely distributed in the global marine environment [1,2,3], which is typically methane gas produced by the escape of gasification from seabed gas hydrates after state instability [4,5]. The seeping methane gas forms bubble plumes that develop in regional aggregations and spill upward from the seafloor, and the escaping gases are mainly greenhouse gases (e.g., CH₄, CO₂, etc.), which may affect global climate change when escaping into the atmosphere [6]. Meanwhile, seafloor gas leakage is also one of the remarkable markers for the identification of gas hydrate formation zones [7].

With the rapid development of ocean exploration technology, the means for seafloor plume detection are becoming more abundant, and the information processing is moving toward multi-source integration. The seafloor gas plume can be observed in an all-round manner by a variety of means such as a sub-bottom profiler, multi-channel seismic exploration, in situ observation, geological sampling, and multi-beam echo sounder system (MBES) [8]. The sub-bottom profiler and multi-channel seismic exploration help to identify anomalous areas of plume seepage by acquiring sedimentary stratigraphic features, but their detection range is limited to a narrow area from the survey vessel route to the seafloor projection, with a low detection efficiency and coverage [9,10]. The new generation of devices that rely on multiple sensors, such as the Innomar SES2000 Quattro or Sixpack, improve the data density of shallow water measurements and can generate 3D palaeohorizons below the seabed. In situ seafloor observation and geological sampling provide crucial data support for the detection of seafloor plumes by directly acquiring acoustic and optical data and geological samples from the target area with high accuracy, but they are limited by the high cost of detection and the difficulty of directly tracing gas seepage points on the seafloor [11]. MBES, characterized by the advantages of a large measurement range, high speed, and high accuracy, has made its traditional version widely available for subsea topographic surveys. With the continuous development of the multibeam system, the backscatter data of MWC and the seafloor can also be collected simultaneously in topography measurement (Figure 1), and MWC detection has the advantages of high coverage, high accuracy, and low cost to become a new choice for seafloor gas plume detection and identification [7,12,13].

The modern MBES is a large integrated system that contains several subsystems with different functions [7,14]. It uses transmitting transducer arrays to send fan-shaped acoustic pulses to the seafloor, and when the sent acoustic waves reach the seafloor through the MWC, receiving transducer arrays receive acoustic pulses and perform beam synthesis. The backscattered acoustic signal is often described by the sonar equation as follows:

$$\mathrm{EL}=\mathrm{SL}-2\mathrm{TL}+\mathrm{TS}$$

(1)

where EL (energy level) is the source level accepted at the multibeam receiver, SL (source level) is the transmitting transducer source level, TL (transmission loss) is the transmission loss accounting for both spreading and absorption, and TS (target strength) is the target strength of the scatterer [7,15]. The obtained water scattering information is able to generate a two-dimensional spatial cross-sectional water column image (WCI) within a single ping [14,16,17]. In the process of multibeam data acquisition, the continuous ping beams launched along the track direction contribute to the construction of a 3D space of the seafloor MWC [13,16] for high-precision and high-resolution 3D observation and analysis of fisheries [18,19], gas emission [20,21], and suspended sediment [22,23] in the MWC. Plume detection, identification, and extraction based on MWC data is currently a hot research topic, and previous authors have mostly analyzed and studied MWC data from the perspective of WCI. The leakage of plume gas in MWC can be detected using WCI, and the methane gas flux at high backscattering flares can be estimated by combining in situ observation and visualization data [20]. Related scholars have achieved plume extraction and recognition within single ping water images from an imaging perspective using various methods such as mask denoising [12], edge detection [24], and deep learning [25]. The images of many factors, such as vessels, suspended matter, signal loss, and sidelobes, during MWC data acquisition lead to the low efficiency and high complexity of gas plume extraction based on image classification and detection angles [26,27].

In the field of surveying and mapping, by means of laser measurement and photogrammetry, the 3D coordinates (XYZ), laser reflection intensity (intensity), and color information (RGB) of each sampled point on the surface of the object under test can be obtained, leading to a collection of points called the “point cloud” [28,29]. The point cloud model is a geometric model with discrete sampling points as primitives, which is a natural representation of the 3D geometric model. The simple data structure and compact storage space of the point cloud model [30] allow for direct representation and processing of optical scan data, and the streamlined and compact data structure allows for the ability to express rich surface detail in 3D geometric models, as no connection topology information needs to be maintained [31,32]. Point cloud datasets are widely used as transitional data models for land and ocean remote sensing [33] and can be combined with geographic information systems (GISs) for analysis, object detection, and 3D modeling [34]. Similar to laser measurements, the bathymetric points obtained from multibeam seafloor topography measurements can be considered as a presentation of 3D point clouds [35,36], and the MBES seafloor topography point cloud data can be used for automatic noise suppression [37] and 3D reconstruction of the seafloor and objects [33]. The high-resolution MWC data obtained by sampling MWC at intervals according to the beam are a record of the acoustic scattering information of the whole MWC from the transducer to the seafloor. By analyzing the data, we can obtain the location, time, reflection intensity, and other information of the sampling point, and it can also be regarded as 3D point cloud data after processing, transformation, and modeling. The 3D geometric information added by the water point cloud model can make up for the sparse sampling points in WCIs and improve the recognition effect of targets in water. However, the MWC point cloud data, which can reflect the morphology of the target in MWC, is also more challenging for the identification and extraction of MWC targets due to its huge amount of sampling points and the presence of more noise.

To extract and identify the seafloor plume, this paper performs spatial homing calculation and other pre-processing for it using the MWC data based on the MWC data features and transforms it into an MWC 3D field point cloud to achieve the point cloud separation and identification of the plume in the MWC, with the feature extraction of the MWC plume based on the 3D point cloud model. Compared with the extraction and recognition methods, from the perspective of imaging, the point cloud model-based extraction of the MWC gas plume can take full advantage of the point cloud data processing to retain more details of the MWC features, more conducive to the rejection of noise signals and the identification of valid MWC targets. The 3D MWC point cloud data model can better help extract the complete geometry of the plume, with higher accuracy and extraction efficiency.

2. Materials and Methods

To achieve the target detection and extraction of the seafloor plume, an overall process is developed as shown in Figure 2. Firstly, the raw multibeam data is parsed to construct a 3D spatial point cloud. Then, the symmetric subtraction and the between-class variance algorithm (Otsu) are employed for MWC noise suppression, and the plume candidate target point cloud is extracted based on the clustering algorithm. Next, the plume feature points are extracted to cluster the candidate target point cloud for feature detection, and finally the target is output to achieve the identification and extraction of the submarine gas plume detected by multibeam.

2.1. Multibeam Data Analysis and MWC Point Cloud Data Location Calculation

MBES for water column imaging has been applied in fishery, physical ocean, and other fields. ELAC Nautik’s Sea Beam deep-water MBES has been widely used, with MWC data used for fish spatial distribution, biomass estimation [38], and bubble detection [20]; Teledyne’s third-generation full bathymetric multibeam system Hydro Sweep DSWCI is equipped with a WCI application—HYDROSTAR WCI Viewer, for online detection of water column targets [39]; and Kongsberg’s EM series deep-water MBES enable a wide range of applications in gas [12,40,41], wreck detection, and fisheries [42]. This article will introduce the multibeam data storage format based on a case study of Kongsberg’s EM series multibeam.

EM series multibeam raw data storage consists of three parts: system parameter module, multibeam data acquisition module, and each external sensor output module. The data are stored in mixed ASCII and binary encoding and split into .all files for bathymetric topographic data, sonar images and system parameters, and .wcd files for MWC data and system parameters. The compatibility module and the file storage module have been added for the storage of its raw data. The following table shows the new EM series high-resolution multibeam datagrams .kmall format processing unit output datagrams (

Table 1. EM datagrams on.kmall format.

Output Datagrams	Datagram Type Code	Description
Installation and runtime datagrams	#IIP	Installation parameters and sensor setup
	#IOP	Runtime parameters as chosen by operator
	#IBE	Built in test (BIST) error report
	#IBR	Built in test (BIST) reply
	#IBS	Built in test (BIST) short reply
Multibeam datagrams	#MRZ	Multibeam (M) raw range (R) and depth (Z) datagram
	#MWC	Multibeam (M) water (W) column (C) datagram
External sensor output datagrams	#SPO	Sensor (S) data for position (PO)
	#SKM	Sensor (S) KM binary sensor format
	#SVP	Sensor (S) data from sound velocity (V) profile (P) or CTD
	#SCL	Sensor (S) data from clock (CL)
	#SDE	Sensor (S) data from depth (DE) sensor
	#SHI	Sensor (S) data for height (HI)
	#SVT	Sensor (S) data for sound velocity (V) at transducer (T)
Compatibility datagrams	#CPO	Compatibility (C) data for position (PO)
	#CHE	Compatibility (C) data for heave (HE)
File datagrams	#FCF	Backscatter calibration (C) file (F) datagram

).

The water backscattering intensity and other data in the MWC datagram are sampled and stored at equal time intervals. It records the backscattering intensity data of a certain area from the seabed to the sea surface with the position of the transducer as the center and the distance from the edge beam to the transducer as the radius. The datagram also records the transmitting angle of each ping data beam, sampling frequency, the number of the sampling point at the beginning of each beam, and the number of the sampling point at the seafloor detected by each beam (0 for absence).

The process of spatial resolution and homing calculation of the MWC sampling point data is shown in Figure 3.

The equation for the spatial homing of MWC sampling points is as follows:

• Calculate the distance between the water sampling point and the transducer:

$$L=\left(idx+Num\right)\times c/\left(2\cdot f\right)$$

(2)

2.

Calculate the water backscattering intensity sampling point along the ship transverse position and water depth:

$$Y=-L\times \cos \left(\mathsf{\alpha }\right)$$

(3)

$$X=-L\times \sin \left(\mathsf{\alpha }\right)$$

(4)

After obtaining the relative coordinates of the MWC backscattered intensity sampling point in steps 1 and 2, the azimuth angle of the beam in the projection coordinates can be calculated, and the absolute coordinates of the MWC backscattered intensity sampling point can be obtained based on the projection coordinates of each ping data transducer.

3.

Calculate the beam azimuth:

$$\mathsf{\beta } =-\mathrm{mod}\left(P\_Head+P\_Con+\mathrm{offset},360\right)$$

(5)

4.

Calculate the projection coordinates of water sampling points:

$$East=P\_East+X\times \cos \mathsf{\alpha }$$

(6)

$$North=P\_North+X\times \sin \mathsf{\alpha }$$

(7)

where $L$ represents the distance between the water sampling point and transducer; idx represents the number of the sampling point of each column receiving beam echo; Num represents the number of the sampling point at the beginning of the receiving beam of each column; $c$ represents the sound speed; $f$ represents the sampling frequency; X and Y represent the transverse position and water depth of the backscattering intensity sampling point along the ship; $\alpha$ represents the beam emission angle; $\beta$ represents the beam azimuth; P_prc represents the projection coordinates (P_ East, P_ North), heading (P_ Head), and meridian convergence Angle (P_ Con) for each ping data acquisition; Offset represents the installation angle of the transducer; (East, North) represents the projection coordinates of water sampling points. Single ping WCI of time-Angle (T-A) and Depth Acrosstrack (D-A) vertical tracks can be generated after analyzing and calculating the data stored in the data packets of external sensors and MWC [42], as shown in Figure 4.

• Time-Angle (T-A) Track Space (Figure 4a): The original point data recorded by MBES. The horizontal axis is marked as the beam number, indicating the beam transmission angle, and the vertical axis represents the sampling sequence. In this display, the bottom is displayed in a near-parabolic shape, and because the beam in the edge direction has its beam footprint significantly widened, the bottom echo signal obtained in that direction is not obvious in its bottom region and presents a widening of the bottom region on the image.
• Depth-Acrosstrack (D-A) Track Space (Figure 4b): The point data after spatial homing. By converting the polar coordinate position of the water sampling point into the absolute coordinate position, the underwater image environment of the water information can be restored more truly.

Cartesian coordinate transformation was carried out for each ping sampling point in the MWC data, and all ping sampling points were reconstructed in space to complete the initial construction of the MWC point cloud space (Figure 4c).

2.2. MWC Noise Suppression

2.2.1. Symmetric Subtraction

MWC point cloud data is huge and contains a large amount of background noise with high values of backward reflection intensity, which is from the diverse sources it contains and is similar to the reflection intensity of the plume [41,43]:

• Noise caused by ships and transducers in MWC will cause an increase in the intensity of acoustic scattering from the surface MWC, which is found in the MWC data within very narrow limits and generally has little impact on the extraction of anomalous targets in MWC.
• A large amount of noise in MWC data is caused by the marine environment, including turbulence, scattering of suspended matters, microbubbles, zooplankton, and other reasons in water, mainly distributed in the strong backscattering layer structure horizontally distributed in the subsurface of seawater, which has a significant influence on the extraction and detection of objects.
• When there is a target with a high scattering intensity or significant changes in seafloor fluctuation, the transmitting beam sidelobe generates noise in WCI with an uncertain location, which is recorded by single or multiple echo sequences in the same time adjacent to each other and presented in the MWC image as a high-brightness arc-like strip. This background noise has strong regularity, which is usually symmetrical with the axis of the central beam [44].

Additionally, the scattering intensity of the sampling points within the MWC data is roughly normally distributed. Since the direct performance of recognition and extraction is computationally intensive, of low efficiency, and prone to recognition errors, noise suppression of the MWC data is required first. According to the symmetric distribution characteristics of noise, the MWC data is demarcated by the central beam for symmetric subtraction after spatial homing to suppress the noise initially for the above noise types, as shown in Figure 5.

2.2.2. Background Point Cloud Removal

Large areas of regular noise within the MWC data return to normal after symmetric subtraction, and a small number of sample point locations with high values of reflection intensity are either targets to be identified or discrete irregular noise points (Figure 6b). At this point, it is necessary to select a reasonable threshold using the Otsu algorithm to segment the background point cloud with a normal reflection intensity of MWC from the point cloud of sites with abnormally high reflection intensity values (Figure 6c).

The Otsu algorithm, also known as the maximum between-cluster variance method, was proposed by the Japanese scholar Otsu in 1979 as an adaptive threshold determination method [45]. The algorithm assumes that the image pixel can be divided into the background and target according to the threshold. Then, the optimal threshold is calculated to distinguish the two types of pixels to maximize the distinction between the two types of pixels.

The calculation is as follows:

• The intensity value of an image is in the range [0, L-1]; then, the threshold k is also in that range so that the threshold k is taken for each intensity value in the range.
• Calculate the between-cluster variance for each threshold.
• Calculate the global optimal segmentation threshold, which is equal to the threshold k that maximizes the variance between clusters; take the mean value in case of more than one threshold.

The results of the Otsu algorithm segmentation are shown in Figure 6.

2.3. Candidate Target Point Cloud Extraction

The segmented field point cloud shows that the plume appears in a high density and small area in space, and the Density-Based Spatial Clustering of Applications with Noise (DBSCAN) algorithm contributes to distinguishing the noise and clustering in the spatial point cloud and extracting them for storage separately [46,47]. The definitions involved in this algorithm are as follows (Figure 7):

• ε-neighborhood (circles): The area within the scan radius (Eps) for a given object.
• Core object (red dot): A point is called a core object if it has no less than the minimum number of contained points (MinPts) in the Eps neighborhood.
• Directly density-reachable (green dot to red dot): For a given dataset D, an object p is said to be directly reachable from object q if point p lies within the Eps neighborhood of point q and q is the core object.
• Density-reachable: A point p is density-reachable from a point q if there is a chain of points p₁, ……, p_i, ……, p_n, p₁ = p, p_n = q, such that each p_{i + 1} is directly density-reachable from p_i.
• Density-connected (yellow dots): A point p is density-connected to a point q if there is a point O such that both p and q are density reachable from O.

Figure 7. Key definition of DBSCAN. (a) Points (green dots) in the ε-neighborhood (circles) are density-reachable from core points (red dot); (b) density-connectivity (yellow dots). Adapted with permission from ref. [48]. Copyright 2012 Copyright Cui.

Figure 7. Key definition of DBSCAN. (a) Points (green dots) in the ε-neighborhood (circles) are density-reachable from core points (red dot); (b) density-connectivity (yellow dots). Adapted with permission from ref. [48]. Copyright 2012 Copyright Cui.

The calculation is as follows:

• Select an appropriate neighborhood value ε and density threshold MinPts according to the characteristics of the data set and mark all data points as unprocessed.
• Randomly select a point P and count the points in the ε-neighborhood of P. If it is greater than or equal to MinPts, mark P as a core point and create a new category. Using P as the starting point, find the points connected to the density of P, and find the maximum set of density connected points. If it is less than MinPts, mark P as a noise point.
• Select another point in the dataset and repeat step 2 until all points are marked as processed.

After the DBSCAN algorithm, the results are shown in Figure 8.

2.4. Plume Target Point Cloud Detection

2.4.1. Model Point Cloud Feature Point Extraction

The actual seafloor plume collected in the National Oceanic and Atmospheric Administration (NOAA) NA080 cruise 2017 in the U.S. Cascadia margin from the open data set of MWC data is used as the model point cloud for this experiment to extract features and learn the plume target detection on the MWC data. Point cloud feature points are those in the point cloud with certain stability and differentiation, and the model point cloud is extracted with feature points for the next step of feature matching. The intrinsic shape signature (ISS) is a feature point extraction method that represents the solid geometric shapes of point clouds. It represents the degree of point identity using the relationship between eigenvalues by calculating the eigenvalue decomposition (EVD) of the scatter matrix of all points in the spherical region [49,50]. The calculation is as follows:

• Create a coordinate system for each point p_i and set the radius r.
• Calculate the Euclidean distance weights for all points p_j in a spherical region:

  $$w_{ij}=\frac{1}{∥p_i-p_j∥},\left|p_i-p_j\right|<r$$

  (8)
• Calculate the covariance array for each point:

  $$\mathrm{cov}\left(p_i\right)=\frac{\sum _{∥p_i-p_j∥<r}{w}_{ij}\left(p_i-p_j\right){\left(p_i-p_j\right)}^{\mathrm{T}}}{\sum _{∥p_i-p_i∥<r}{w}_{ij}}$$

  (9)
• Calculate the eigenvalues $\left\{{\lambda }_i^1,{\lambda }_i^2,{\lambda }_i^3\right\}$ of each point $\mathrm{cov}\left(p_i\right)$ and rank the eigenvalues from largest to smallest.
• Set the threshold ε₁ and ε₂, and the point satisfying the following conditions is the ISS feature point:

  $$\frac{{\lambda }_i^2}{{\lambda }_i^1}⩽{\epsilon }_1,\frac{{\lambda }_i^3}{{\lambda }_i^2}⩽{\epsilon }_2$$

  (10)

  where p_i represents the point picked by the algorithm, p_j represents all points in the spherical region calculated by p_i, r represents the spherical region search radius, $w_{ij}$represents the Euclidean distance of all points in the spherical region, $\left\{{\lambda }_i^1,{\lambda }_i^2,{\lambda }_i^3\right\}$are the eigenvalues of the points, and ε₁ and ε₂ are the thresholds of the ISS feature points.

2.4.2. Target Point Cloud Recognition Based on FPFH Features

The overall geometric characterization of the point cloud can be achieved using the point cloud feature points for their neighborhood feature calculation. The fast point feature histogram (FPFH) is a simplified form of the point feature histogram (PFH) computational approach, which is to parameterize the spatial differences between query points and neighboring points and form a multidimensional histogram to describe the geometric properties in the k-neighborhood of a point (Figure 9a). Instead of computing the combination of all adjacent points, FPFH simplifies its computational complexity by computing the simplified point feature histogram (SPFH) for each point in the k-neighborhood of the query point, and then weighting all the SPFHs into the final fast point feature histogram [51,52] (Figure 9b).

The calculation is shown as follows:

• By establishing a local coordinate system to calculate the relative relationship between two points p_t and p_s in the k-neighborhood, a local coordinate system with u, v, and w as axes is created with n_t and n_s as the corresponding normal vectors, and p_s as the coordinate origin. The calculation is as follows:

$$\{\begin{array}{l}\mathrm{u}={\mathrm{n}}_{\mathrm{s}}\hfill \\ \mathrm{v}=\mathrm{u}\times \frac{\left({\mathrm{p}}_{\mathrm{t}}-{\mathrm{p}}_{\mathrm{s}}\right)}{∥{\mathrm{p}}_{\mathrm{t}}-{\mathrm{p}}_{\mathrm{s}}∥{}_2}\hfill \\ \mathrm{w}=\mathrm{u}\times \mathrm{v}\hfill \end{array}$$

(11)

2.

Compute the angular variations of n_t and n_s as follows:

$$\begin{array}{ll}\alpha \hfill & =\mathrm{v}\cdot {\mathrm{n}}_{\mathrm{t}}\hfill \\ \varphi \hfill & =\mathrm{u}\cdot \frac{\left({\mathrm{p}}_{\mathrm{t}}-{\mathrm{p}}_{\mathrm{s}}\right)}{d}\hfill \\ \theta \hfill & =\arctan \left(\mathrm{w}\cdot {\mathrm{n}}_{\mathrm{t}},\mathrm{u}\cdot {\mathrm{n}}_{\mathrm{t}}\right)\hfill \end{array}$$

(12)

where D represents the Euclidean distance between two points. By calculating the $\left(a,\varphi ,\theta ,d\right)$ quaternions for each point pair in the neighborhood, it is possible to reduce the coordinate information of the two points and their normal vectors associated with 12 parameters (xyz coordinate values and the corresponding normal vector components) to four. This results in a relative relationship between two points that can be used to describe surface features.

3.

At the end of the calculation, the quadratic parameters between all point pairs are put into the histogram in a statistical way. Each parameter is divided into b subintervals, and finally, the number of subintervals falling into each subinterval is counted separately to obtain the PFH feature descriptors.

FPFH for the query point only computes a set of $\left(\alpha ,\varphi ,\theta \right)$ between that point and the points in the neighborhood to obtain SPFH.

The formula is as follows:

$$FPFH(p)=SPFH(p)+\frac{1}{k}{\sum }_{i=1}^k\left(\frac{1}{w_k}·SPFH\left(p_k\right)\right)$$

(13)

where w_k represents the weight and the FPFH feature histogram stores feature vectors of 33 elements. Figure 10 shows the visualized FPFH feature histogram of two feature points a and b in the selected model point cloud.

FPFH is used to identify the matching features between point clouds, and the key points of the model and candidate point clouds are detected for features by setting the spatial relation threshold and matching threshold to determine whether they are target point clouds [53]. The thresholds defined for feature matching are as follows:

The matching threshold consists of scalars in the range [0, 1]. Two feature vectors are matched when the normalized Euclidean distance is less than or equal to the matching threshold. The spatial relation threshold, on the other hand, uses point cloud data to estimate the spatial relationship between points associated with potential feature matches, and finds potential matches based on the spatial relationship threshold.

For FPFH feature matching, the normalized Euclidean distance matching feature is calculated based on the found matching points, the algorithm initialization detection threshold is set to 0.1, and the scores are returned as a column vector to evaluate the degree of matching, where a smaller value indicates a higher probability that the point cloud to be detected is a plume. The threshold is adjusted to the optimal matching state in the experiment and determines the threshold, as shown in Figure 11. The scores are averaged to determine whether the cloud at the detection point is a plume.

3. Results

3.1. Data Acquisition

In this paper, three experimental areas were selected for method validation using multibeam data collected by the Kongsberg EM302 multibeam system aboard the NOAA Okeanos Explorer in the northern Gulf of Mexico in late summer 2011 (Figure 12). The cruise carried a 30 kHz Kongsberg EM302 multibeam system, a 18 kHz Kongsberg EK 60 single-beam system, and collected backscattered intensity from the MWC and seafloor in the northern Gulf of Mexico.

The ship used the onboard Applanix POS/MV (ver. 4) to record and correct multibeam data for any ship’s motion before being logged by SIS software. The C-NAV GPS satellite service system provided DGPS correctors to the POS/MV with the positional accuracy expected to be better than 2.0 m. All the corrections (motion, sound speed profile, sound speed at sonar head, draft, sensor offsets) were applied during real-time data acquisition in Kongsberg data acquisition software Seafloor Information System (SIS) ver. 3.8.3. Sippican XBT casts (Deep Blue, max depth 760 m) were taken every 6 h.

3.2. Process and Results of Experimental

By analyzing the MWC data in the experimental area for noise suppression, we extracted a total of some candidate target point clouds as shown in Figure 13, which were irregular noises within the plume or MWC.

Then, the target detection of the candidate point clouds was performed through the method and threshold selection in Section 2.4 of this paper, candidate point clouds with scores less than 0.1 were selected, and finally, the plume target detection results were output, verifying the location of the plume seepage nozzle compared with those observed in real time during the cruise (Figure 14).

3.3. Experimental Conclusion

Experimental area 1 is located at the edge of the salt dunes in the Gulf of Mexico, and 3D seismic data in this area indicate its past carbonate or hydrate formations, justifying the large-scale gas leakage in experimental area 1. By comparing the voyage measured gas leakage points and the acoustic profiler observation study of the water plume in this area, the salt dune edge in the northern Gulf of Mexico is the active area of methane plume leakage, and other areas also have sporadic plume leakage points [54,55], which are consistent with the detection results of this method, and the practicability of the method is verified.

The proposed method relies on the construction of point cloud space to achieve MWC noise suppression and plume recognition and extraction. The two plumes extracted by the detection of the method in this paper are shown in Figure 15. The construction of the point cloud space preserves the partial structure of the plume, which cannot be determined in the single ping image, and more fully preserves the spatial structure of the plume.

Compared with the single ping identification and extraction method, the recognition efficiency is enhanced, and the recognition process is simplified. The test time of three experimental areas is shown in the following

Table 2. The time of the algorithm process.

	Length of Route Line (m)	Number of Ping	Time of Noise Suppression (s)	Time of Point Cloud Extraction (s)	Time of Target Detection (s)
area 1	17,441.15	1130	39.29	4.56	10.13
area 2	2469.47	178	6.56	3.16	7.17
area 3	16,197.88	1402	30.98	1.84	8.23

. The length of the route line, ping number, and time of each step of the experiment were counted in the three test areas. The noise suppression in experimental area 1 cost 39.29 s, which was the longest step in each step. The time of point cloud extraction and identification is relatively random, which is related to the noise and the number of plumes in the experimental area. If there is less noise in MWC, it only takes a short time to extract and identify the plume point clouds, as shown in the experimental area 3.

4. Discussion

This method builds a 3D point cloud scene model of MWC based on the acoustic reflection intensity, extracts and detects objects in MWC through the 3D geometric features of the point cloud, and verifies the feasibility of the method through experimental areas. In the process of target detection and extraction using MWC data, the accuracy of the method is affected by the data itself and many problems in the process of constructing the point cloud space. its performance may be influenced by the following factors.

4.1. Characteristics of MWC Noise

In this paper, using the methods of symmetric subtraction and the Otsu algorithm, we mainly suppress the ocean environment and the sidelobe interference noise and extract the anomalous point cloud in MWC. As shown in Figure 16, the range of the sound scattering intensity of MWC is about (−120 dB, 8 dB) before symmetric subtraction, with a peak value of −60 dB, compared to about (−70 dB, 75 dB) after symmetric subtraction, with a peak value of 0 dB. The data volume at the peak increases from $3.7 \times$10⁶ to $4.0 \times$10⁶. After the Otsu algorithm, the amount of sound scattering intensity data in MWC is greatly reduced, in a double-peak structure, with the peak values located at −24 and 10 dB, respectively. According to the statistics of the plume scattering intensity data volume subsequently extracted, the peak of 10 dB is at the high probability of the plume distribution while the peak of −24 B is the anomalous noise within MWC.

The histogram of the water reflection intensity at different stages of the proposed method allows for more intuitive observation and analysis of the characteristic relationship between the MWC noise and plume reflection intensity (

Table 3. Distribution patterns of the backscattered intensity in MWC at different stages of the algorithm.

	Original Reflection Intensity	After Symmetric Subtraction’s Reflection Intensity	After Otsu’s Reflection Intensity	Gas Plume Reflection Intensity
Data volume	186,662,016	186,662,016	9301	2612
Median	−38	0	3	11
Skewness	−0.3329	0	−0.4827	−0.5836
Kurtosis	3.2696	3.6584	3.1081	3.9414

). The original water reflection intensity showed a negative skewness distribution with a peak value of 3.2696. After symmetric subtraction, the water reflection intensity showed a standard normal distribution and its peak value increased to 3.6584. It can be seen that the method of symmetric subtraction in this paper has a significant effect, and the large-scale abnormal noise points were mostly caused by the marine environment. After the Otsu algorithm, the amount of data decreased dramatically, and the noise points with high reflective intensity that were not removed temporarily in the water led to the double-peak structure presented in the histogram. According to the mean value, the skewness and peakness of the backscattering intensity of the final plume reached the maximum value. It can also be seen that the extracted plume noise points are basically removed in this paper, and the plume itself has a highly reflective intensity and a negatively skewed and spiked distribution.

4.2. Spatial Distribution of MWC Point Cloud

In this paper, based on the reflection intensity information of the sampling points of MWC data, spatial point cloud modeling was carried out to realize the plume identification and extraction of the 3D point cloud. In the extraction of 3D point clouds, the hierarchical clustering algorithm based on the Euclidean distance uses KD-tree to search the neighborhood and classifies point clouds according to the spatial Euclidean distance between proximal points as a judgment criterion. The calculation is shown as follows:

• Select the seed points; the KD-tree is used to search the radius R neighborhood of the seed points. If there are points in the neighborhood, they are grouped into the same cluster Q.
• Select a new seed point in cluster Q and continue to perform step (1). If the number of points in Q does not increase, Q clustering ends
• Set the threshold interval of cluster points. If the number of points in cluster Q is within the threshold interval, the clustering results will be saved.
• Select a new seed point from the remaining point cloud Q and continue to perform the above steps until all points in the point cloud are traversed.

The above clustering methods are compared and analyzed with the results of the DBSCAN algorithm in this paper. As shown in Figure 17, it is found that the hierarchical clustering algorithm has a good clustering effect, but it can only classify the point cloud by the cluster, and its segmentation effect is not satisfactory for the small discrete noise points around a certain cluster. In contrast, the DBSCAN algorithm has a better clustering effect. It can independently identify MWC noise in the point cloud based on the density to remove it, and it can more accurately identify discrete points around the point cloud cluster as noise removal. The comparison results show that consideration of the water target and water density difference between the noise points of the DBSCAN algorithm is unaffected by the clustering points set K. The arbitrary-shaped high-density points can find the water glowed class, in terms of both plume flow water extraction and anomaly target discovery in water extract, and also has great application value.

The construction of the MWC 3D space relies on the stacking of per-ping data by the track direction, using the acoustic reflection intensity information to respond to the target object within MWC [12]. However, the resolution of MWC data according to the track direction is much lower than that of sampling points within ping, which leads to faultiness and sparsity in the water space, and has adverse effects on the target extraction, subsequent gas flux estimation, and 3D reconstruction of the target. Currently related scholars use Fledermaus, Voxler 3, and other software to visualize the target objects in the MWC data in a virtual 3D environment with volume rendering and 3D data display [56] and construct WCI as a 3D voxel grid to reduce the effect of MWC data discontinuity and improve the accuracy of subsequent MWC target processing [12,40]. Improving the spatial inhomogeneity of MWC point clouds will improve the accuracy of target identification and detection of MWC data and provide a strong basis for subsequent work.

4.3. Experimental Results outside MSR

When the transmitted beam reaches the seabed, its reflection intensity is much higher than the backscattering intensity, and the reflected wave reaches the receiving transducer first. Due to the presence of the receiving beam sidelobe, the echo is recorded by most beams, thus forming a strong interference with the same time, which is reflected in the MWC data as a half-arc-shaped strong reflection zone. This reflection zone is called specular reflection, and the radius of this arc is called the minimum slant range (MSR) [41]. The MWC data within MSR are all from seawater and of high quality while the data outside MSR contain significant interference from the receiving sidelobe [12]. The algorithm used in this paper achieved a good recognition effect for the data within MSR, but the application of the algorithm outside MSR has some limitations. In the noise suppression stage, the seabed part outside the MSR range cannot be completely removed, which has a significant impact on the subsequent point cloud extraction, as shown in Figure 18.

Although using the data inside the MSR for target detection can improve the efficiency and accuracy of target recognition, it significantly weakens the advantage of the large coverage width of MBES and limits the application of MWC data. Therefore, it is of great significance to improve the quality of data outside MSR and improve the utilization rate of MWC data.

5. Conclusions

This paper proposes a method to extract and identify the gas plume based on the MWC point cloud by calculating the spatial homing of MWC data, transforming the MWC data into a 3D field point cloud, and realizing 3D imaging, feature extraction, and identification of the seafloor plume based on the spatial point cloud. According to the analysis of MWC noise, the MWC noise suppression method was proposed to realize the initial segmentation of abnormal targets in the MWC point cloud. According to the spatial distribution characteristics of the MWC point cloud model, the DBSCAN algorithm and FPFH feature detection algorithm were proposed to realize the extraction and detection of 3D targets of the MWC plume. Compared with existing methods, this study reduced the difficulty of plume contour shape identification, improved the performance and accuracy of plume identification and extraction based on the 3D point cloud, and verified the reliability of this method using data from the northern Gulf of Mexico.