Technical Note : Detection of gas bubble leakage via correlation of water column multibeam images

Introduction Conclusions References


Introduction
Sonar offers great potential for remote sensing of the oceans because underwater transmissions of acoustic signals are much more efficient than optical signals.Aside from seafloor mapping, echosounders are established as standard tools for the remote detection of targets in the water column, such as submarines, single fish, fish shoals, zooplankton, or gas bubbles.Small wind-induced or ship-entrained surficial gas bubbles have been acoustically investigated in detail (Leighton, 1994;Medwin and Clay, 1998).Early publications focused on the remote acoustic sensing of methane gas emissions released from the seafloor into the water column (McCartney and Bary, 1965;Merewether et al., 1985).Such methane gas emissions were reported to occur world-wide from virtually all continental margins, estuaries, and river deltas (Judd and Hovland, 2007).The magnitude of these seafloor methane gas emissions and the consideration of potential bubblemediated flux to the atmosphere, where methane acts as a very strong greenhouse gas, is not well understood (Kvenvolden and Rogers, 2005).The importance of gas bubble sensing in the future is underscored by the need to monitor and evaluate potential gas leakage associated with carbon capture and storage (CCS) projects (Oldenburg and Lewicki, 2006).
Progress in digital signal processing brought about the design of so called multibeam systems (de Moustier, 1988), covering large angles (160 • ) and delivering high resolution (0.5 • beam angle).Schneider von Deimling et al. (2007) and Nikolovska et al. (2008) demonstrated its potential for gas seepage mapping.Modern water column imaging (WCI) multibeam systems are used for obstacle avoidance, bioacoustic investigations, coastal navigation assurance and for Published by Copernicus Publications on behalf of the European Geosciences Union.scientific seafloor mapping (Gerlotto et al., 2000;Mayer et al., 2002).Best et al. (2008) first applied particle imaging velocimetry (PIV) cross-correlation techniques to multibeam water column data in order to extract alluvial flow structures.
The use of multibeam data for WCI is hindered by the inconvenient handling of very large amounts of data.Effective algorithms are required for low false-alarm target detection and target classification.The cross-correlation approach presented in this paper aims to extract unique data patterns that can only be produced by rising gas bubbles.Such patterns were derived from an example dataset acquired in the Baltic Sea at 24 m water depth using a next generation WCI 50 kHz multibeam system.The proposed processing scheme will be useful to design a very sensitive, reliable, and automated seep/leak bubble detector.

Gas bubble acoustics
Gas bubbles in water generally very efficient scatter acoustic energy.The reasons for this are the very large differences in sound speed and density between water and gas, as well as resonance effects controlled by bubble size, frequency, and water depth (Minnaert, 1933;Medwin and Clay, 1998).In our field test (24 m water depth, 50 kHz frequency) we expected scattering in the geometric regime with target strengths (TS) between −65.7 and −40.1 dB, calculated using the non-modal approximation (Medwin and Clay, 1998) for spherical gas bubbles between 1 mm and 20 mm diameter.Thus, even the smallest gas bubbles are expected to cause sufficient backscatter for detection.
For more information on gas bubble hydroacoustics, the reader is referred to Medwin and Clay (1998), Leighton (1994), andAnderson (1950), who give the full solution to sound scattering from a fluid sphere in liquid by solving the three dimensional wave equation.

Multibeam
An L-3 ELAC Nautik prototype SB3050 was installed below the moonpool of R/V Poseidon with small, but mobile, 3 • × 2 • transducers in a Mills-Cross configuration.The SB3050 system comes with the sophisticated WCI software HYDROSTAR WCIViewer, that was used for online water column inspection.The system transmits three frequencycoded sectors that cover a swath of 140 • with 315 equiangle beams.Full water column imaging is available almost in real-time, including geo-referenced amplitude backscatter values for each beam b and respective samples s (7 cm).The system was operated with 0.15 ms, 50 kHz continuous wave pulses and a ping repetition rate of 0.15 seconds.Ship motion (roll, pitch, heave, yaw) and ray refraction was compensated to obtain depth z and athwart distance x relative to the transducer.Additionally, the system performed "roll stabilization", giving rise to a receive-swath that is directed perpendicular to the seafloor (center beam = incidence angle) for any roll motion up to ±10 • .Range-dependent corrections, accounting for acoustic absorption and geometrical spreading, were achieved by application of a time-variant gain.To cover the respective backscatter intensity for gas bubble visualization we selected an optimized color scale and video recording of the data (no raw data access).Finally, sonar bitmap images (640 × 433 pixels, 8-bit) with pixel indices x and z were generated (MB(x , z ) t ) for each ping time t.

Automated bubble detection via PIV
For the proposed bubble detection algorithm, PIV was integrated into the post-processing workflow for the analysis of multibeam sonar images.This method applies a crosscorrelation of intensity values of images taken at a given time interval to generate velocity vectors of individual "seeds" within the images.Bubble-induced backscatter anomalies are one such seed type in sonar images.For more details of PIV methods the reader is referred to, for example, Raffel et al. (1998).
In Fig. 1, the schematic workflow of the proposed bubble detection algorithm is illustrated.Two sonar bitmap images MB(x ,z ) t1 and MB(x ,z ) t2, separated by the ping interval t (in our data t = 0.15 s), are compared via crosscorrelation.
Each bitmap is subdivided into a 9 by 6 matrix (indices i,j ) where each element represents a corresponding subwindow or "interrogation window" (Fig. 1a).For every interrogation window, the intensity values at t 1 are compared to the adjacent values at t 2 , within a certain displacement ( x, z; Fig. 1b).Maximum correlation values are then selected (Fig. 1c) to determine the prevailing velocity vector in the respective interrogation window i,j , which is most influenced by moving gas bubbles, if present.Since the seeds (bubbles) were not evenly distributed in the data, a normalized crosscorrelation was applied.
For a better performance, cross-correlation was calculated in the frequency domain (achieved with the MATLAB mpiv toolbox), where the cross-correlation was transformed into a complex conjugate multiplication, thus reducing the computational effort from O [N 4  ] to O[N 2 log 2 N] operations.Prior to this, inter-frame mean values were subtracted in order to suppress noise with noise-noise or signal-noise mismatch.Subsequent inverse Fourier transformation gives the 2-Dlag-correlogram (Fig. 1c) containing the normalized correlation coefficient R for a given displacement x and z within the time interval t.Peak finding (Fig. 1c) was performed with a Gaussian peak-finding algorithm implemented in mpiv.Incorporation of water current and bubble rise velocity estimates allow for the definition of an "approved zone" of expected bubble lag displacements (Fig. 1b, c; black rectangle) as an accept/reject parameter.Finally, the gathered pixel lag displacements x and z were converted into distances (meters) x and z by the pixel-to-distance relation for the x and z directions.Division by the ping interval t results in velocity vectors v i,j (x,z) in m/s for each interrogation window (i,j ).

Field Site
In the Western Baltic Sea, vast areas are characterized by shallow gas located 2-4 m below the seafloor in the "Holocene mud" (Laier and Jensen, 2007).Hydroacoustic evidence of very shallow gas in this layer was presented in Fiedler and Wever (1997) at the study site.Here, a rosette-water sampler system (RWS) was deployed onto the Holocene mud to cause some overpressure on to the gascharged sediments at 24 m water depth.Subsequently, rising gas bubbles in the water column were imaged for 75 s, at a fixed position (54 • 15.109 N, 11 • 31.942E).

Time slice
In Fig. 2, the backscatter values MB(x, z) t1−3 between incidence angles of −43 • and +44 • at three different times are presented as color coded "polar" images.Low backscattering plotted in blue shades, whereas high backscattering is plotted yellow-red shades.High backscatter caused by the seafloor is manifested in the prominent semi-circular side-lobe echoes across all beams.Another artifact persistently plots in the center beam with elevated backscatter values caused by simultaneous noise affecting the center beam of the planar array -mostly by constructive interference after beamforming.
Figure 2a shows the downcast of the RWS in the form of a strong backscatter signal in the water column, accompanied by side-lobe-echoes in all beams.
A few seconds after the gear had landed individual moving patterns emerge (Fig. 2) which are interpreted to represent rising gas bubbles squeezed out by the load of the RWS on the gas-charged sediments.The bubbles were slightly shifted to the right during the tracking period of a few seconds between Fig. 2b and c.Approximately one minute after deployment the bubble pattern disappeared from the sonar images.

Beam slice
In our experiment, the vessel stayed at a fixed position and therefore subsequent pings span a 3 vector space temporally.To investigate the rising velocities of gas bubbles, we need to analyze their temporal evolution.This is achieved by resorting the 3-D data cube of consecutive pings in single beam slice images.Each beam slice represents a single beam, now imaging travel time against subsequent pings (Fig. 3, beam 31 with 500 pings), where rising gas bubbles emerge very clearly as straight lines.
During the downcast of the RWS (ping 0-180), some horizontal lines appear between pings 170-280 (Fig. 3 I).This pattern is barely visible in Fig. 2 and is regarded as neutrallybuoyant gas bubbles detaching from the down-going RWS.Subsequently, rising straight lines plot between pings 330-500 and bubble rise velocity can be determined from z/ t (corrected for projecting the travel time of Beam 31 to the center beam/vertical direction) in Fig. 3   the steepness of the gradient and the cease in backscatter intensity indicates smaller bubbles with lower rise velocity and reduced backscattering strength TS.
This kind of data presentation highlights the bubble trajectories in each beam (Fig. 3), but also improves visual object detection, because of the signal coherency in consecutive pings.If lateral currents and bubble drift are low compared to their rise velocity, then a fixed beam over the entire sampling period should be preferred for the center beams.

Correlation Processing
To sense moving targets in the sonar images, PIV was applied to determine the dominating velocities therein.The resulting velocity vectors are presented in Fig. 4. Interrogation windows (5, 4) and (5, 5) show distinct negative (upward) v z .In (7, 3) and (9, 4) some artifacts show up at the boundary between the swath image and the background image.However, they do not occur consecutively.Outside of the bubbleinfluenced area (off-seep), only low magnitude vectors appear in the respective interrogation windows To investigate data over a longer period, v z and v h values were calculated between successive pairs of images (p = 300 to p = 400) and visualized in a compound vector plot for an "on-seepage" interrogation window (5, 4) and an "offseepage" (Fig. 5) interrogation window (6, 3).The (5, 4) "on-seepage" vectors are strongly biased towards +x and −z directions, in agreement with bubble movement observed in Fig. 2 and Fig. 3.The corresponding lateral displacement values of +0.2 m s −1 could not be confirmed because current information is missing.
The off-seepage analysis reveals neither a bias towards +x/−z direction, nor to smaller values of v z within the seepindicative v z seep range (defined later)."Off-seep" vectors plotted in Fig. 5 are thus interpreted as derived from noisenoise correlation.Green and red rectangles indicate the "on" and "off" seepage interrogation window, respectively.
Fig. 5. Velocity vector components v z and v h calculated for the "on-seep" interrogation window (5, 4) and "off-seep" interrogation window (6, 3) in green and red, respectively.For the location of the interrogation window, see Fig. 4.

Seep bubble detection
Natural gas seepage can be discriminated into minor and major gas release (Leifer and Boles, 2005).Minor gas release is characterized by individual or groups of gas bubbles following each other with a narrow Gaussian-distributed bubble size spectrum.Typical peak diameters were reported at ∼5 mm at depths up to 1000 m (Leifer and McDonald, 2003;Leifer and Boles, 2005;Sahling et al., 2009;Sauter et al., 2006;Schneider v. Deimling et al., 2011).Terminal rise velocities of such buoyant gas bubbles can be estimated after Clift et al. (1978), e.g. 1 mm to 20 mm large bubbles would rise at a v z seep rate between 0.1 and 0.35 m s −1 .In contrast, major gas releases are driven by overpressure and cause broad bubble spectra ranging between 0.3 mm-50 mm and upwelling flows yield v z values exceeding 2 m s −1 (Leifer, 2006).
To our knowledge, gas and/or oil seepage/leakage from the seafloor represents the only mechanism to sustain a continuous production of rising objects with v z seep in the ocean.Thus, our detection scheme aims to find characteristic "rise patterns" in the data, which are regarded as seep-indicators (v z seep ).Fish represent a main source of confusion with respect to acoustic bubble investigations, but these confusions can be overcome by consideration of differences in echo trajectories of bubbles and fish (Ostrovsky, 2003 and2009).Although fish may rise temporarily within the v z seep range, persistent detection of v z seep values clearly points to a gas source on the seafloor.
In particular, oceanographic current data could improve the proposed processing scheme (Fig. 1) through velocity investigations on passively rising bubbles, as pointed out in Schneider von Deimling et al. (2010).If current data are not available, a best assumption about the local velocities could be made.
Stationary acoustic gas bubble monitoring studies (Dworski and Jackson, 1994) have concentrated on elevated backscatter from bubbles -not their trajectories.However pure backscatter analysis as a gas bubble seep indicator is very vulnerable to fish echoes.
Maximum and mean rise velocities of v z min =−38 cm s −1 and v z =−16 cm s −1 were resolved with our PIV processing.Both of these values fit very well to the manually picked rise velocities (Fig. 3II, -34 cm s −1 ) and theoretical values for minor seep bubble spectra.However, the 0.15 ms pulse limits the vertical resolution to approximately 0.1 m, and thus the ping to ping resolution is small.At least for low flux gas leakages, the rising behavior of gas bubbles lies in a limited range that can be used for discrimination between rising gas bubbles and other unwanted echo targets.It must be stated here that our study represents a proof of concept with limited resolution.Much more accuracy will be possible with a next generation high frequency system, e.g.operating between 400 kHz and 2 MHz (supplemental video).
Improved sonar imaging or signal processing is needed to detect gas bubble leakage.Difficulties arise with increasing beam angle, where the prominent semi-circular first arrival side-lobe (Fig. 2a) strongly disturbs the water column signal in modern WCI multibeam mapping sonar data.Here, sophisticated methods such as "optical flow" may be valuable for improved detection of moving patterns.In comparison to image processing, signal processing of raw WCI data is very www.ocean-sci.net/8/175/2012/Ocean Sci., 8, 175-181, 2012 computationally expensive.Therefore, a great deal of pioneering work will be needed to handle the large quantity of WCI data.Nevertheless, basic signal processing of raw WCI data allowed us to geometrically resort our test data.Our data were recorded while the vessel was stationary, and bubbles were not significantly affected by lateral currents perpendicular to the swath.However, the trajectories of echoes caused by buoyant gas bubbles in the water column are three dimensional in nature.Therefore, a single swath can only insonify a two dimensional slice (snapshot) of this three dimensional volume.By collecting multiple swaths over a period of time, a data volume of subsequent pings can be stored.After resorting these data to horizontal time/depth slices, a vertically contiguous object is cut into a series of horizontal slices, in which it is then imaged as an isolated object.Thus, it is possible to trace the spatially-coherent bubble rise pattern down to the seafloor, even if currents have deflected the gas bubbles from one swath into another.
Future multibeam mapping systems may offer athwart swath steering of about 20 • .In combination with a high ping rate, this would allow for repeated insonification of individual gas bubbles during a survey.
Methods are already available to extract the third component of the velocity vector as well (stereo techniques, dual plane PIV: Raffel et al., 1998), and thus allow for tracking a 3-D bubble path.However, the computational requirements for 3-D-PIV become even higher.

Limitations
One major constraint in using the proposed correlation processing arises from seeding density and resolution issues.Both the horizontal and vertical resolution are strongly dependent on the range and frequency.If the distance between individual or groups of bubbles (seeds) is significantly smaller than the data resolution, then the acoustic trace of gas seepage plots spatially continuously and the introduced interrogation windows are no longer seeded with individual particles.In this case, PIV might fail.This would be particularly relevant in case of a major gas release with variously sized bubbles having very different rise/dissolution behavior.
The prototype system in our study is not ideally suited for such tasks given its relatively low horizontal resolution (2 • × 3 • , 50 kHz) and large pulse length (0.15 ms) which limits vertical rise measurements.Today, multibeam mapping systems and imaging sonars are commonly available operating at frequencies up to, and even beyond, 400 kHz.Application of such high resolution systems would offer millimeter to centimeter vertical resolution at high ping rates.Under these conditions, the proposed processing scheme is expected to work much better and a quantitative assessment would become feasible e.g. for short range applications such as AUV surveys or in situ monitoring.
A promising step would be the direct application of the detection algorithm on the signal data.Although in this study we have only used signal processing for re-sorting the data to generate beam-slice images, future developments and algorithms could be applied to the raw data directly.

Conclusions
With the advent of the next generation WCI multibeam systems, water column investigations become realistic and require new processing techniques.With data of our 50 kHz prototype sounder we successfully developed a method to record the rise pattern of gas bubbles released from the seafloor.The applied method is based on correlation techniques adapted from image processing and potentially represents a basis for automated seepage/leakage detection at short ranges.Correlation processing of raw WCI data, in this study tested and only used for geometrical aspects, is a promising tool in future WCI investigations, as the methods adapted from image processing can be directly applied to raw data.With higher frequency systems and sector steering techniques the method will be more sensitive.Our study shows that automated WCI processing in general can potentially be a major asset for a wide range of hydroacoustic applications, such as gas seepage analysis and monitoring of man-made subsea installations, with special emphasis on low gas flux leakage detection and ocean current analysis.
Acknowledgements.The authors would like to thank L3-Communications ELAC-Nautik GmbH for supplying the prototype echosounder as well as for system adaption, installation, and support within the framework of the German SUGAR project (grant no 03G819A).Special thanks go to the scientific diver team namely Florian Huber, Christian Howe, Nikolaus Bigalke and Björn Thoma for tricky transducer installation and to the crew members of R/V Poseidon for their support.We honor the review of two anonymous reviewers and the critical proof reading of Ralf Prien and Michael Glockzin.The ship cruise and respective research leading to these results has received funding from the European Community's Seventh Framework Program (FP/2007(FP/ -2013) ) under grant agreement n • 217246 (Baltic Gas).Supplemental video data could be acquired by the support of R2Sonic, Ageotec, and Embient GmbH supplying a R2Sonic 2024 with WCI under the framework of ECO2 (grant no 265847).

Fig. 1 .
Fig. 1.Schematic diagram illustrating the workflow of a proposed detection algorithm showing (a) synthetic echo-images with a solitary "seed" (black rectangle) moving into +x and −z direction over time period t, (b) cross-correlation in the frequency domain and the definition of the "approved zone" controlled by horizontal currents v x and rise velocity v z , and (c) bubble displacement visible as a Dirac delta peak within the "approved zone" in a 2-D correlogram.

Fig. 2 .
Fig. 2. Color-coded raw backscatter images recorded by the multibeam echo sounder showing background water column reflections (dark blue), center beam artifacts and semi-circular side-lobe-interference close to the seafloor (light blue).(a) The downcast of a rosette-watersampler (yellowish RWS) is seen.(b) Gas bubble release (yellowish ellipses).(c) Vertically shifted gas bubbles.The swath width x(z) at depth covers 41 m, where z represents the water depth minus transducer sea surface installation offset (3.6 m).

Fig. 3 .
Fig. 3. "Beam slice" representation with the x axis corresponding to ping times (1 ping = 150 ms) and the y-axis showing two-waytravel time for Beam 31.The downcast of the rosette-water-sampler (RWS) is seen in the left of the image.Horizontal straight lines are seen at p = 150 (I).Further to the right of the image (at p ∼ 320), a compound of rising gas bubbles occurs (II).

Fig. 4 .
Fig. 4. Mean subtracted image similar to Fig. 2 superimposed with the pixel displacement vector (green arrow) composed of x and z displacements calculated for each interrogation window (i,j ).Green and red rectangles indicate the "on" and "off" seepage interrogation window, respectively.