Pseudo-3D cubes from densely spaced subbottom profiles via projection onto convex sets interpolation: An open-source workflow applied to a pockmark field Open Access

Subsurface geologic features are 3D in nature and thus often require a tailored 3D imaging approach to enable reliable interpretations (Yilmaz, 2001). Such 3D seismic reflection surveys using relatively low frequencies (approximately 1–500 Hz) are widely used in hydrocarbon exploration; however, due to practical and financial constraints, 2D seismic surveys are still much more common in academic research. There is a growing interest in high-resolution 3D imaging of the shallow seafloor subsurface, especially for offshore wind farms (e.g., Samuel, 2022; PGS, 2023), which requires modifications to the acquisition and processing methods used in hydrocarbon exploration. Several dedicated 3D systems for high-resolution applications have been developed since the late 1990s, e.g., the P-Cable system (Planke et al., 2009). These systems consist of high-resolution seismic sources (e.g., air gun, boomer, and sparker) with towed streamers (Müller et al., 2002; Monrigal et al., 2017), whereas other systems use semirigid hydrophone arrays (SEAMAP-3D) (Müller et al., 2009, 2013). Alternative systems for imaging the shallow subsurface are subbottom profilers (SBPs), often using chirp or sweep frequency sources. Their benefits are wide bandwidths, high frequencies (kilohertz range), and controlled-source signatures that facilitate high vertical resolutions (centimeter-decimeter range), an improved signal-to-noise ratio (S/N), and quantitative analyses (Quinn et al., 1998; Gutowski et al., 2002). Such systems use transducer and receiver arrays that are either hull mounted (Kongsberg Geoacoustics Ltd., 2019) or arranged within a rigid array configuration (Gutowski et al., 2002; Bull et al., 2005).

Two-dimensional SBP data are commonly recorded during conventional lower-resolution academic and industry voyages as either simultaneously collected ancillary or distinct data sets for a specific reason. Their main applications include archaeological studies, geotechnical site investigations, and shallow geologic mapping (e.g., Bull et al., 2005; Müller et al., 2009, 2013). For many research targets or areas, however, dedicated 3D acquisition systems are often not available or applicable. For seismic applications, 3D imaging can be achieved by exploiting the effect of strong surface currents leading to streamer feathering (e.g., Lin et al., 2019) or using model-guided interpolation for regional data sets (e.g., Trinchero et al., 2014; O’Keefe et al., 2017; Kim et al., 2019). Alternatively, pseudo-3D volumes of the subsurface have been generated by combining dense grids of 2D profiles followed by 2.5D or 3D poststack migration (e.g., Zhao et al., 2023). In contrast, very few examples of 3D imaging using SBPs exist. For example, Kim et al. (2020) use a conventional chirp SBP towfish mounted to the ship’s hull to detect and image small buried targets in very shallow water depths (<10 m) by creating a pseudo-3D volume from 39 parallel 2D profiles. However, for deeper water beyond coastal areas, only a single study that we are aware of describes the application of a hull-mounted chirp SBP system to generate a pseudo-3D volume from 60 profiles (Kim et al., 2016). There, a flex binning approach was used to infill empty bins by incorporating traces from neighboring bins (Kim et al., 2016).

During the research voyage TAN2006 with R/V Tangaroa in 2020 (Pecher et al., 2022), densely spaced high-resolution SBP lines were acquired using the hull-mounted Kongsberg Maritime TOPAS PS18 parametric SBP (Kongsberg Geoacoustics Ltd., 2019) across a paleo-pockmark field on the Chatham Rise off the east coast of New Zealand’s South Island (Figure 1). Here, we present a novel open-source workflow to generate comprehensive pseudo-3D cubes imaging subsurface features in high resolution. In contrast to the approach of Kim et al. (2016), we use an interpolation method to infill empty bins after sparse 3D binning resulting from the difference in average line spacing (approximately 25 m) and trace distance along the line (approximately 5 m). We aim to preserve the highest spatial resolution possible and enable subsequent analyses of geologic features in the shallow subsurface, which is required for the previously mentioned applications such as wind farm site surveys.

In marine reflection seismic acquisition, the recorded data are often irregularly sampled in the spatial domain due to physical or economic conditions such as varying ship speed, weather conditions, time constraints, or other effects such as streamer feathering and dead channels. However, nonequispaced and undersampled data can severely affect subsequent processing steps such as surface-related multiple attenuation (e.g., Gao et al., 2010, 2013) or wave equation-based migration (Wang et al., 2014). Therefore, regularizing and interpolating seismic data is an essential processing step (e.g., Wang et al., 2015a; Wang, 2016; Jicheng et al., 2018).

In general, 3D seismic volumes can be subdivided into data sets that are (1) regularly sampled, (2) irregularly sampled, or (3) sparse (Chen et al., 2018). Regularly sampled describes data on an equidistant grid with periodically missing traces, whereas purely irregular data in the spatial domain can be transformed into a regular grid structure with randomly missing traces. Reconstructing highly sparse data sets poses a significantly greater challenge (Chen et al., 2018) and generally requires additional information, for example, to compensate aliasing (Ghaderpour, 2019). At present, methods for irregular seismic data interpolation can be classified into several categories: wave equation-based (e.g., Spitz, 1991; Trad et al., 2002; Naghizadeh and Sacchi, 2007; Liu and Fomel, 2011), signal processing-based (e.g., Xu et al., 2005; Ghaderpour et al., 2018; Jicheng et al., 2018; Dai et al., 2021; Cao et al., 2022), rank reduction-based (e.g., Trickett et al., 2010; Oropeza and Sacchi, 2011; Ma, 2013; Jia et al., 2016; Bayati and Trad, 2023), geostatistical (e.g., Turco et al., 2019; Yeeh et al., 2020), and deep-learning methods (e.g., Jia and Ma, 2017; Wang et al., 2020).

Signal processing-based methods, such as the antileakage Fourier transform (Xu et al., 2005), fast iterative adaptive approach (Dai et al., 2021), antileakage least-squares spectral analysis (ALLSSA) (Ghaderpour et al., 2018) and its adaption, the multichannel ALLSSA (Ghaderpour, 2019), or projection onto convex sets (POCS) (Abma and Kabir, 2006) and its modifications, have been widely used and improved over the past few decades due to their efficiency and overall simplicity (Zhang et al., 2017). Among the most popular approaches are sparsity-promoting methods closely related to the compressed sensing paradigm (Candès et al., 2006; Dai et al., 2021). They follow the general assumption that seismic signals are sparse in a transform domain (Dai et al., 2021; Zhao et al., 2021). Data interpolation is formulated as an inversion problem (Cao et al., 2022), which can be solved using least-squares minimization solutions and convex sparse inversion models (e.g., Yang et al., 2013; Ghaderpour et al., 2018; Cao et al., 2022). One benefit is that no a priori information is required as the interpolation is solely based on the seismic signal (Zhao et al., 2021). Many single-scale or multi-scale mathematical transforms capable of reconstructing sparse data have been introduced during the past few decades. Examples of single-scale (i.e., global) transforms are the fast Fourier transform (FFT) (Xu et al., 2005; Abma and Kabir, 2006; Gao et al., 2010, 2013; Ghaderpour et al., 2018; Ghaderpour, 2019), Radon transform (Trad et al., 2002), and wavelet transform (Zhao et al., 2021). Multi-scale transforms include the dreamlet transform (Wang et al., 2014, 2015b, 2015a), seislet transform (Fomel and Liu, 2010; Gan et al., 2015, 2016), shearlet transform (Häuser and Ma, 2012; Jicheng et al., 2018), and curvelet transform (Candès et al., 2006; Yang et al., 2012; Cao and Zhao, 2015; Zhang et al., 2015, 2017, 2020). Single-scale transforms are generally less effective at characterizing seismic data; however, they are computationally more efficient. Multi-scale transforms exhibit better sparse expression potential (Zhao et al., 2021) due to their multidirectionality and anisotropy, allowing for the representation of 2D curve-like features (e.g., Candès et al., 2006; Yang et al., 2012; Cao and Zhao, 2015).

Among the available algorithms to solve the optimization problem, the POCS algorithm is generally preferred due to its simplicity, robustness, computational efficiency, and interpolation quality. Conceptually, it is a two-step iterative method, which obtains the estimated solution in the transform domain before reinserting the observed data to improve the convergence (Wang et al., 2015b; Wang, 2016; Wang and Lu, 2017). The POCS algorithm was first introduced by Bregman (1965) and used to interpolate missing traces in 3D seismic data (Abma and Kabir, 2006) using a simple linear, iteration-dependent threshold model. Subsequent studies established exponential (Gao et al., 2010; Ge et al., 2015; Zhang et al., 2015, 2017; Zhao et al., 2021) and data-driven threshold models (Gao et al., 2013), which improved the convergence and, thus, the computational efficiency of the POCS algorithm. An adaptive method to calculate the minimum threshold based on the sparse coefficients was recently proposed by Zhao et al. (2021). Several studies indicated an increase in efficiency and interpolation quality by transforming the input data from the time to the frequency domain prior to the actual sparse transform (Zhang et al., 2015; Wang, 2016; Wang and Lu, 2017; Zhao et al., 2021). Wang (2016) exploits the symmetry of real-valued input signals in the frequency domain using only positive frequencies for interpolation.

One deficit of the original POCS algorithm is the inability to simultaneously interpolate data and remove random noise from the sampled traces, as the noisy original data are reinserted at each iteration step (Abma and Kabir, 2006). Gao et al. (2013) propose a weighted POCS (WPOCS) to reduce random noise using a weighting factor $α∈(0,1]$⁠, which was subsequently improved by exchanging the order of the iterative threshold operator and the projection (i.e., reinsertion) operator (Wang et al., 2015b; Wang, 2016). The fast POCS (FPOCS) was proposed by Yang et al. (2013) and Gan et al. (2015, 2016) using data from a previous iteration step to precondition the following one, showing significant performance improvements due to a much faster convergence. For efficiency reasons, almost all 3D implementations of POCS variants use a time or frequency slice approach for trace interpolation (Abma and Kabir, 2006; Cao and Zhao, 2015; Zhang et al., 2015; Kim et al., 2018), which is also implemented in our workflow (Figure 2).

The study area is located on the western Chatham Rise off New Zealand’s South Island in water depths between 550 and 580 m (Figure 1). The Chatham Rise is a submerged component of the Zealandia continent (Mortimer et al., 2017) and is framed by the deep Hikurangi Plateau to the north and the Bounty Trough to the south (Figure 1b). It is a relatively shallow (approximately 400–700 m water depth) flat-topped ridge, extending >1000 km eastward from the South Island. Large parts of the Chatham Rise seafloor are covered with pockmarks of varying sizes (Davy et al., 2010; Stott et al., 2019), which seem to be bathymetrically controlled and, for the study area, restricted to water depths of approximately 450–700 m (Stott et al., 2019; Pecher et al., 2022). For the study area, Davy et al. (2010) describe abundant “textbook” pockmarks on the modern seafloor (Figure 1a) and paleo-pockmarks buried in the subsurface.

Prior to the 3D interpolation of the recorded TOPAS profiles, several preprocessing steps are required to ensure consistency, remove noise bursts, and, more importantly, correct potential misties at intersections of the 2D profiles. Due to the high vertical resolution of SBP systems, the latter effect is substantially more pronounced than for multichannel seismic data with lower frequencies. To generate a comprehensive 3D cube from the dense TOPAS grid, we developed a novel open-source processing workflow consisting of two stages: (1) the (pre)processing of each 2D profile and (2) the generation and subsequent interpolation of a 3D cube (Figure 3). All processing steps are implemented using the programming language Python and the well-established packages “NumPy” (Harris et al., 2020), “pandas” (The Pandas Development Team, 2020), “SciPy” (Virtanen et al., 2020), “segyio” (Equinor ASA and The Developer Team, 2017), “Xarray” (Hoyer and Hamman, 2017), “SEGY-SAK” (SEGY-SAK Developer Team, 2020), and “dask” (Dask Development Team, 2016), among others. Please refer to Appendix  B for a list of the Python packages used here, also available as requirement files (see the “Data and materials availability” section).

This study uses more than 1575 km of SBP data, arranged as a dense grid of >200 lines with approximately 25 m spacing (Figure 1a). This data set was recorded using the Kongsberg Maritime TOPAS PS18 SBP (Pecher et al., 2022), a hull-mounted nonlinear parametric echosounder using the simultaneous transmission of two slightly different high-frequency signals to generate difference frequencies (Kongsberg Geoacoustics Ltd., 2019). The benefits of this method are a wide bandwidth and a narrow beam with no sidelobes, resulting in very high vertical resolutions (<15 cm) while achieving relatively high subsurface penetration (up to 200 m). During TAN2006, a linear frequency-modulated chirp signal of 20 ms length with a return bandwidth of 4 kHz (2–6 kHz) was triggered in conjunction with the simultaneously operating multibeam echosounder (MBES) EM302. The transmitted beam is motion stabilized and accounts for the seafloor slope using multibeam depth information. To reduce recorded file sizes while retaining a high sampling interval (here 0.025 ms), SBP data are commonly recorded in vertical time windows. Window lengths of either 200 or 250 ms were used during TAN2006, which resulted in 10,000 and 12,500 total samples per trace, respectively, given the selected sampling interval (Pecher et al., 2022). The applied sound speed velocity was inferred from a hull-mounted probe and typically varied from 1495 to 1505 m/s during operation. Using an equispaced operation mode and a survey speed of approximately 2.6 m/s (5 kn), the average trace spacing is approximately 5 m for most profiles with a beam footprint of approximately 45 m for 560 m water depth (4.5° beamwidth) in the study area. In addition to a receiver gain, an initial high-pass filter (>1 kHz), and a matched filter, no additional processing steps were applied before the recorded data were stored in SEG-Y format (Pecher et al., 2022).

Small SEG-Y files (<50 traces), a consequence of line name changes during acquisition, were merged with the previously or later recorded files, depending on the spatial distance and the consecutive field recording number in the trace headers. The trace positions during acquisition were stored as geographical coordinates (WGS84, EPSG:4326) but transformed into the projected coordinate reference system WGS84/UTM zone 60S (EPSG:32760) to facilitate geometric operations during subsequent processing steps.

When changing the delay recording time (delrt) of the specified recording window while recording an SEG-Y file, an incorrect delrt was occasionally assigned to the last and first traces before or after the adjustment, respectively. These delrt needed to be corrected before compensating for different delay times within a file. A Python script was developed to (1) detect changes in delrt based on trace header information and then (2) compare neighboring traces using the position of their maximum amplitudes within a defined window to identify incorrect times (Figure 3, delrt correction). Subsequently, the corrected SEG-Y files must be compensated for different delay recording times to remove significant vertical offsets within one file (Figure 3, trace padding). The minimum and maximum recording times of all traces were extracted before each trace was padded with zeros at the start or end, respectively. Without the previous correction of incorrect delay times, such traces would appear as outliers at incorrect vertical positions.

Another environmental effect to consider is the varying tidal elevation during the acquisition, which extended over several days (Figure 3, tide compensation). The observed sinusoidal tidal amplitude in the study area is approximately ±0.8 m or ±1 ms (Figure 5a), resulting in noticeable vertical offsets, especially for neighboring profiles recorded at different times on the tidal cycle. As the profile lengths (5–10 km) and thus acquisition times (<1 h) are relatively short, the tidal variation along a single profile shows a linear trend (Figure 5b). Consequently, the tidal compensation detrends the profiles resulting in a continuously increasing or decreasing deviation from the original depth along each recorded profile (Figure 5c). To compensate for this effect, we developed a stand-alone Python package called “tpxo-tide-prediction” (Warnke, 2021) that provides a minimal interface to predict tidal elevations based on the TPXO-9 atlas models (Egbert and Erofeeva, 2002). For each TOPAS profile, the time-dependent tidal elevation with reference to the mean sea level is computed, and each trace is statically adjusted accordingly.

Even after applying the previous processing steps, vertical offsets might persist at profile intersections, possibly caused by insufficient static correction. In general, misties in seismo-acoustic data might originate from different amplitudes, phase rotations, or time shifts (Bishop and Nunns, 1994). In this study, we only considered the latter effect because all profiles were acquired using the same system and configuration (Figure 3, mistie correction). According to Harper (1991), the profiles in this study represent a closed-loop configuration, where a least-squares approach can be used to solve the set of overdetermined equations (Ghaderpour et al., 2018). Our implemented Python algorithm (1) detects all intersections (Figure 6b), (2) loads k nearest traces from associated profiles and calculates average reference traces (Figure 6c and 6d), (3) extracts a user-defined window and calculates the trace envelope (Figure 6e), (4) computes the crosscorrelation of both reference traces (Figure 6f), (5) filters low-quality (i.e., noisy) intersections by setting a threshold based on the Pearson’s correlation coefficient, and (6) computes the global (i.e., line-wise) mistie (i.e., time shift) by minimizing the residuals via the least-squares solution.

All TOPAS profiles show a varying intensity of distinctive noise bursts of approximately 20 ms. Interestingly, their abundance varies during the survey, and their signal characteristics indicate an interference of subsequently emitted pings with the recorded trace. We developed a Python algorithm to remove most noise bursts (Figure 3, 2D despiking). Without this despiking, such spikes would significantly impact the quality of later interpolation, as they introduce artificially high amplitudes at random locations in space and time. A user-defined moving window approach is used, where the center traces are compared with the neighboring ones. These windows may overlap by a certain percentage on the time axis (here 25%), whereas the step size in the spatial domain is set to one trace to ensure that each trace is filtered. Spikes are defined as samples where the reference (center) trace amplitude exceeds the chosen mean, median, or root-mean-square (rms) background amplitude (as per the user’s choice). All flagged spike samples are merged and potential misdetections (e.g., a small number of consecutive samples) are removed. The flagged samples are then either scaled to the background amplitude level or replaced with zeros, or the mean, minimum, maximum, or median amplitudes; the latter were used for this study. The quality of the despiking algorithm largely depends on selecting a suitable window size and the empirically determined threshold parameter. Here, we use a window of 30 ms by five traces and a threshold factor of four, a configuration that effectively removes most of the noise bursts.

For better input/output (I/O) performance and to enable later parallel processing using “dask”, the preprocessed TOPAS profiles are converted from SEG-Y to the widely adopted netCDF-4 format (Rew and Davis, 1990) using the Python package “SEGY-SAK” (SEGY-SAK Developer Team, 2020). All relevant metadata information and setup parameters for initializing the sparse output netCDF volume are stored in human-readable YAML configuration files. In the first step, the parameter files are loaded and a temporary geometry is created, where the inline orientation is aligned in a north–south direction. If the rectangle sides defined by the input corner points are not aligned accordingly, the corner points are then rotated by a specific angle (in degree) around the rotation center using an affine transformation that preserves the collinearity of the input (Figure 7a). Next, the trace locations from all 2D profiles are loaded, either from previously created auxiliary files or directly read from the SEG-Y headers. The trace locations, output geometry, and bin sizes are then used to define the inline (referred to as “iline”) and crossline (“xline”) bin ranges of the output before each trace is assigned to an appropriate inline/xline bin (Figure 7b). For single bins containing multiple traces, we tested four approaches: computing (1) the mean or (2) median amplitudes, (3) selecting the trace that is spatially nearest to the bin center, or (4) applying an inverse distance weighting (IDW) algorithm to create a reference trace (Figure 7c). The IDW algorithm was selected to produce the sparse cubes presented in this study.

Gain recovery is essential to homogenize traces from different origins and compensate for amplitude decay with distance. The amplitude of elastic waves emitted during seismo-acoustic surveys is affected by several factors when propagating through a medium. Among these, spherical divergence and absorption have the most significant impact (Yilmaz, 2001). Spherical divergence, the decrease of amplitude A with distance due to geometric spreading, can be compensated for by applying a simple time-variant gain function $A(t)gain=A(t)*t2$ (Figure 3). This approach does not correct for the absorption but preserves the relative amplitudes of the seismic data. We used a custom Python implementation of the gain calculations defined in the “sugain” module of the Seismic Unix software package (Cohen and Stockwell, 2022) to compensate for the spherical divergence. Traces from different data sets might also show varying amplitude intensities due to differing acquisition settings, which is especially common for multichannel seismic data. Although all data in this study have been acquired with the same configuration parameters, environmental conditions such as the sea state might influence the recorded signal strength. Therefore, we balanced all traces in the sparse cube with reference to their respective rms amplitude $A(t)b=A(t)/rms$ (Figure 3, amplitude balancing).

The TOPAS data were recorded with a sampling interval (⁠$Δt$⁠) of 0.025 ms. Because the bandwidth of the return chirp signal is limited to 2–6 kHz (Figure 8a) and a frequency spectrum analysis shows that the main frequency content is 2–5 kHz, the data can be considered to be oversampled. To reduce processing time in the subsequent interpolation, we downsample each trace by a factor of two resulting in a new sampling interval (⁠$Δtr$⁠) of 0.05 ms without affecting the main signal frequency content (Figure 8b). The original Nyquist frequency (⁠$fN$⁠) of 20 kHz was reduced to 10 kHz through downsampling, which is about twice the maximum signal frequency of approximately 5 kHz (Figure 8b). Prior to resampling, the data were band-pass filtered to remove low- and high-frequency noise outside the signal bandwidth (corner frequencies: 1800/1900/5000/5200 Hz; Figure 8, lower panels). SBP data sets are commonly displayed using the phase-independent trace envelope (i.e., instantaneous amplitude) (Henkart, 2006) because of its better interpretability compared to the correlated waveform that often shows a more “ringy” appearance (Quinn et al., 1998; Saustrup et al., 2019). Especially for legacy data sets, often only the envelope data are available. Therefore, we calculated the trace envelope using the Hilbert transform, shifting the signal to a lower frequency range (<2 kHz; Figure 8c) to show the applicability of the proposed workflow for more widespread data types.

The core part of our workflow is the interpolation of sparse 3D data that can be conducted in the time or frequency domain (Figure 2). However, the computational efficiency of the interpolation is significantly improved when converting the input data into the frequency domain first by applying the FFT along the time axis (see the “Projection onto convex sets” section and Appendix  A). Retaining only positive frequencies for interpolation effectively reduces the computational load by a factor of two (Wang, 2016; Ghaderpour et al., 2018). This data reduction does not affect the result because the input signal consists of only positive real 32 bit floating point values (float32). We applied an additional low-pass filter (5000/5300 Hz) based on the frequency spectrum analysis after the domain conversion (Figure 8c, lower panel) to further improve the interpolation runtime. Frequencies above the low-pass filter cutoff are excluded from the interpolation. During each iteration, the POCS algorithm applies a 2D forward transform to each frequency (or time) slice before executing hard or soft thresholding (Figure 2). The updated sparse coefficients are input for the subsequent inverse transform. In the last step of each iteration, the observed data $Dobs$ are reinserted (Figure 2). If the maximum number of iterations is reached or the cost function error is below a user-defined threshold, the iteration is stopped (Figure 2). This slice-wise interpolation of sparse input data poses an embarrassingly parallel problem as the processing of each frequency (or time) slice is entirely independent of adjacent ones. To accelerate the computational efficiency, we use the open-source framework “dask” (Dask Development Team, 2016) to parallelize the computation by interpolating multiple slices simultaneously using different central processing units (CPUs). This approach further allows scaling the workload from a single machine to distributed clusters without significantly changing the code basis. Each process outputs a batch of n interpolated slices as a netCDF file before all batches are merged into a single file. Finally, the inverse FFT is computed along the frequency axis to transform the reconstructed data back into the time domain after applying an optional Gaussian smoothing filter over the interpolated slices to reduce residual noise, otherwise resulting in a less sophisticated appearance (Figure 2).

As explained in the “Pseudo-3D processing workflow” section, one of the, if not the most, crucial prerequisites for generating a pseudo-3D cube is to account for vertical offsets at line intersections. A static correction is primarily applied to improve the reflection coherence along a single 2D profile, but the vertical repositioning of the trace in time also affects observed offsets at intersections. Although the absolute shift is relatively small (±0.2 ms; Figure 9a), applying a static correction significantly increases the output image quality (Figure 10). However, when plotting the detected seafloor reflector for each trace, distinct vertical offsets are visible for individual profiles (e.g., black arrows in Figure 9b). These offsets are less pronounced after compensating the data set for the changing tidal elevation ranging between approximately ±1 ms (Figure 9c). Comparing the seafloor depth before (Figure 9b) and after the tide compensation (Figure 9d), fewer profiles with an apparent vertical offset are observable in the data, although some profiles still display a noticeable downward shift (the black arrows in, Figure 9) and depict large relative misties (the dark red lines in Figure 9e). After successful mistie correction and, thus, removal of previously observed artifacts, the resulting seafloor appears much smoother and homogeneous (Figure 9b and 9d). The apparent residual offsets in the western map sections result from the chosen visualization method (the colored points with transparency) as well as varying line spacing (see Figure 1) and not real uncorrected offsets.

A comparison between a raw unprocessed TOPAS profile and its fully processed equivalent is shown in Figure 10. The unprocessed data clearly show abundant noise bursts of approximately 20 ms in length (the left inset plots, Figure 10a and 10b) and a sawtooth-like shape of the seafloor reflection (the right inset plots, Figure 10a and 10b). After execution of all (pre)processing steps, most noise bursts have been removed, and the static effect impacting the reflection coherence has been corrected. Consequently, the output profiles appear less noisy and fragmented, thus facilitating better interpretability (Figure 10c and 10d). The full waveform (Figure 10a and 10c) and the later interpolated trace envelope data sets (Figure 10b and 10d) are visualized for comparison.

A critical factor affecting the interpolation quality is the selection of the most effective sparse transform. Although multi-scale transforms (e.g., shearlet and curvelet) should generally result in better interpolation quality than single-scale transforms (e.g., FFT or wavelet transforms), this is at the expense of computational efficiency, as indicated in Table 1. A visual comparison of interpolated slices in the time and frequency domains for four different sparse transforms depicts evident quality variations (Figure 11), supported by the sharpness index (Table 1). For this comparison, we used the default POCS parameters (Table 2) but set $pmin$ to a constant value of 0.0001 and selected the sparse transform accordingly. Adjusting the remaining parameter only slightly altered the results (Figure 11), which are representative of the following observations. The FFT shows promising results with some random noise in time and frequency slices (Figure 11a and 11b), also visible in the wiggle trace representation (Figure 11i and 11j, black lines). In contrast, the wavelet transform using a Coiflet with a scaling function coefficient c = 5 (“coif5”) fails to reconstruct the missing information for several bins (Figure 11c and 11d), resulting in artificial high or low reconstructed amplitudes (Figure 11i and 11j, the blue lines). The sharpness index of the wavelet transform is higher than for the FFT, but the apparent holes render the result unusable. Similar outputs have been generated for different wavelet types that show similarly poor quality and are therefore not visualized here.

Surprisingly, the interpolation quality of the multi-scale transforms is lower than that of the FFT (Figure 11e–11h), also reflected in the calculated sharpness indices (Table 1). The shearlet-based interpolation exhibits a prominent pattern of intertwined linear artifacts covering the time slice, resulting in a blurred appearance (Figure 11e) and incorrect amplitudes (Figure 11i and 11j, the green lines). Moreover, it fails to reconstruct the missing data in the frequency domain as clearly visible artifacts in a grid alignment of higher relative amplitude can be identified (Figure 11f). Linear artifacts are also identifiable but less pronounced using the curvelet transform (Figure 11g). The effect is reduced in the frequency domain but still noticeable along the slice boundaries (Figure 11h). One benefit of the curvelet transform over the FFT is its ability to effectively reduce the noise content in the output, also shown for the reconstructed traces (Figure 11i and 11j, the red lines). This transform is provided by the CurveLab package (Candès et al., 2006); however, it is only available for Unix-like operating systems. In addition, the runtime for the multi-scale curvelet transform is considerably longer than the single-scale FFT while producing comparable results.

Using the FFT as the sparse transform, we generated interpolated pseudo-3D cubes in the time and the frequency domains (Figure 12) using a bin size of 5 m × 5 m and the default POCS parameters (Table 2). When comparing both cubes, it is evident that the time-domain interpolations show less pronounced features in the time slice (Figure 12b) than in the frequency-domain interpolation (Figure 12c). The power spectra of interpolated slices at 775 ms two-way traveltime (TWT) and sharpness indices (Table 3) indicate a slightly better interpolation result for the frequency domain. At the same time, vertical and horizontal linear features of relatively higher amplitude are more evident for the time domain (Figure 12b). These artifacts are caused by the grid cell orientation and the low data coverage of the input data prior to interpolation (approximately 32%), as apparent from the sparse cube (Figure 12a). In addition, the upper and lower sections of the selected inline and crossline profiles feature more noise with relatively high amplitude in the time domain (Figure 12b). However, the most significant advantage of the frequency-domain interpolation is its overall runtime of 8 min 14 s. In contrast, the runtime was 50 min 11 s for the time domain, which is surprisingly longer than the projected execution time inferred from single-slice experiments (Table 1). A more in-depth analysis showed that the interpolation took substantially more time on average (>10 s) compared with approximately 1.12 s per slice used for estimations (Table 1). Some of this can be explained by the computational overhead from file I/O, but the runtime is still about twice as long as for the frequency domain (approximately 5.6 s). In addition to the fact that only half the frequency slices need to be interpolated, the total number of slices to process is further reduced after applying a low-pass filter, as described in the “Stage 2: Sparse data interpolation (3D)” section. Consequently, computational efficiency and interpolation quality observations favored the use of the frequency-domain approach, which is in agreement with previously published results (e.g., Wang, 2016; Wang and Lu, 2017; Zhao et al., 2021).

In addition, numerous parameter tests revealed that adjusting the total number of iterations is one of the most impactful variables in determining the interpolation quality of the output cube with a fixed bin size (e.g., 5 m × 5 m). Using only 10 iterations results in vertical stripes in the profile sections and a distinct cross pattern in the power spectrum of time slices, representing the grid cell boundaries (data not shown). Small improvements are visible when using 50 or 100 total iterations (data not shown), with the sharpness indices indicating minor improvements for increasing iterations (Table 3). Because an iteration-dependent thresholding model is used as part of the POCS algorithm, 50 iterations are sufficient to generate relatively stable intermediate results with only minor adjustments of the sparse coefficients between subsequent iteration steps (Figure 13a). To visualize this effect, the metric defined in equation A-11 is plotted for each iteration step of the four experiments (Figure 13a). For a better comparison of these experiments, three frequencies are investigated as examples. Although the absolute error declines with the increasing number of iterations, the interpolation of high frequencies depicts a random behavior (2500 and 5000 Hz, Figure 13a). In contrast, the 15 Hz slice shows an exponential decay as expected when using an exponential threshold model (Table 2). These patterns are expected as the spatial variability of lower frequencies is much less compared with higher frequencies, which explains the described variability for the first half of the total iterations indicated by the error functions (Figure 13a). Furthermore, the computational runtimes scale almost linearly with increasing iterations (Figure 13b). This chart also emphasizes the benefits of a parallel implementation as a serial execution of, e.g., 50 total iterations would increase the computation time by almost a factor of approximately 16, which is equivalent to the number of threads (16) used for the parallel interpolation (Figure 13b).

An appropriate bin size largely depends on the acquisition setting and, more specifically, on the original trace spacing along and the line space between available input 2D profiles. All presented cubes were binned using an output grid cell size of 5 m × 5 m (xline × iline), which is roughly equivalent to the average trace spacing and naturally results in a sparse 3D cube considering the line spacing of approximately 25 m (Figures 1 and 12a). The initial coverage for a bin size of 5 m × 5 m is approximately 32%, thus leaving approximately 2/3 of the bins empty and requiring interpolation. Increasing the bin size by factors of 1.5 (7.5 m × 7.5 m) and 2 (10 m × 10 m) not only leads to better bin coverages and a higher sharpness index (Table 3) but also removes some details from the interpolated results (Figure 12a and 12b). If the bin size is increased to 10 m × 10 m, the residual effect of the bin cell boundaries is noticeably reduced when compared to a higher resolution binning, where these artifacts are still slightly visible in the power spectrum (Figures 13c and 14a). For an even larger bin size of 15 m × 15 m (Figure 14c), the interpolated cube shows a blurred appearance where finer structures and also noise are suppressed. An alternative option to using a sophisticated interpolation workflow would be binning using a line spacing of 25 m and filling the remaining gaps using a straight-forward nearest neighbor or, as applied in this example, linear interpolation (Figure 14d). Despite the larger bin size, this approach preserves more detail in the generated output but at the cost of reflection coherence. Furthermore, noise is introduced to the binned traces, resulting in a fragmented appearance (Figure 14d). This observation is supported when examining the power spectrum of the noisy time slice that shows high amplitudes scattered throughout the spectrum, indicative of changes at high spatial frequencies due to abundant noise.

The average line spacing of the selected example area within the study region is reasonably consistent; however, it is about twice as large for the profiles to the west and east (Figure 1). Therefore, we assessed the capability of our workflow to deal with a predominant acquisition direction without perpendicular lines, a common pattern for SBP data recorded during 3D seismic surveys using, e.g., the P-Cable system (Planke et al., 2009). We interpolated the sparse 3D cube using only profiles in north–south and west–east orientations, filling primarily the inlines and crosslines of the sparse cube, respectively (Figure 15b and 15c). Of course, this artificial data reduction results in significantly lower bin coverages of approximately 20.7% and 14.3% for only the inlines and crosslines, respectively (Figure 15d). The interpolated cube using only inlines shows smearing along the predominant north–south direction in the presented time slice and vertical artifacts in the example crossline 600 (Figure 15b). This directionality is also indicated in the power spectrum of the times slice, which exhibits a relatively higher amplitude for the normalized spatial frequencies for the inlines than the crosslines (Figure 15b). The equivalent observation can be made for the data reduction using only crosslines for interpolation (Figure 15c). There, the vertical artifacts in the inline sections are more pronounced because of the lower initial bin coverage. This results in a fragmented appearance, where many of the more small-scale features, such as seafloor pockmarks, are not resolved (xline: 800, Figure 15a and 15c).

One possibility to account for the unequal trace spacing within and between the lines is to increase the bin size perpendicular to the predominant profile orientation (Figure 15e and 15f). If only inlines are available, adjusting the bin size in the inline direction by a factor of three to 5 × 15 (xline × iline) reduces the number of empty bins between the north–south-orientated profiles. This results in an initial bin coverage of approximately 57%, which is in agreement with the expected threefold increase from the original coverage (Figure 15d). Naturally, the displayed cube depicts a lower spatial resolution in the crossline direction; however, the previously visible smearing in the time slice and fragmentation in the profile section are suppressed (Figure 15d), resulting in better image quality (Table 3). Comparable observations can be derived upon examination of the interpolated cube using a larger bin size in the inline direction when limiting the input data set to only crosslines (approximately 41% coverage) (Figure 15d and 15f). Because individual traces are assigned a higher weight with a decreasing number of bins, this might result in preserving noisy traces, as shown between xlines 225 and 250 (Figure 15f), which are only indistinctly recognizable when using a homogeneous bin size (5 m × 5 m, Figure 15a).

Residual vertical offsets and misties in 2D profiles pose a critical constraint on data quality following 3D interpolation. We used the STA/LTA algorithm for a custom seafloor reflection detection followed by a data-based statistical approach to identify and remove misdetections before computing the residual static for each trace. In contrast to our profile-based method, Kim et al. (2016) apply a five-point smoothing filter on an extracted horizon and use the calculated offsets to remove the static effect. However, this method requires a complete, not sparse, 3D cube to pick the reference horizon manually and is, therefore, unsuitable for our interpolation-based workflow. Another option is to use MBES bathymetry as a reference surface to tie in the detected seafloor samples (Kim et al., 2017, 2020). As part of this workflow, we provide a generic workflow for the static correction solely based on the input 2D SBP data. In areas showing localized and abrupt morphological changes such as the observed seafloor pockmarks (Figures 1, 4, 5, 9, and 10), satisfactory results were produced by limiting the maximum vertical shift of the affected traces after manual inspection (Figures 4 and 10).

Unlike many previous studies, the multi-scale curvelet (e.g., Yang et al., 2012; Cao and Zhao, 2015; Zhang et al., 2015, 2017, 2020; Wang and Lu, 2017) and shearlet transforms (e.g., Häuser and Ma, 2012; Jicheng et al., 2018) proved to be not only less computationally efficient (Table 1) but also achieved lower interpolation quality (Figure 11; Table 1). These results are a little surprising because the global FFT is capable of reconstructing linear and quasilinear events but requires a windowing approach to locally approximate curved events (Gao et al., 2010). One possible explanation might be that our algorithm interpolates individual frequency (or time) slices rather than seismic shot gathers or 2D poststack sections in the time-space (t-x) domain, which often feature curved and hyperbolic events. In contrast, especially the frequency slices display more fuzzy patches of high amplitudes rather than the distinct quasilinear and curvilinear reflectors commonly known from seismic sections (Figure 11). Compared with the local curvelet and shearlet transforms, the FFT proved to be more robust against relatively noisy and sparse input data. In general, numerical experiments on synthetic and 2D field data indicate an improved S/N for interpolations using multi-scale transforms, particularly when performing the interpolation in the frequency domain (Wang, 2016; Jicheng et al., 2018; Zhao et al., 2021). In our case, the minor addition of random noise by using the FFT is preferable over the insufficient reconstruction of the multi-scale transforms (Figure 11). A general problem with implementing an efficient POCS algorithm is that most published studies focus on synthetic and relatively simple 2D examples and only occasionally 3D field examples, but all of these comprise low-frequency seismic data. As one of the only published studies that show interpolation of sparse 3D (land) data, Kim et al. (2018) demonstrate satisfactory results with a curvelet-based POCS. In their study, they implement an intermediate step before applying the curvelet transform by transforming the input slice into the $kx−ky$ domain, which enhances the output quality but almost triples the runtime (Kim et al., 2018). Based on our extensive parameter tests, the performance loss and the additional layer of dependency when using shearlet (Lee, 2016) and curvelet transforms (da Costa, 2021) in Python, we opt for the FFT as the most robust and effective transform, which is also part of de facto standard scientific Python packages such as “NumPy” (Harris et al., 2020) and “SciPy” (Virtanen et al., 2020).

As the numerical experiments have shown, using a higher number of total iterations per slice is typically desirable to achieve a higher interpolation quality (Figure 13). However, the improvements are negligible above a certain number of iterations. For this study, we used 50 iterations as a compromise between interpolation quality and runtime (Figure 13). Fewer iterations led to more persistent noise, whereas more iterations increased the runtime without noticeable quality improvements. The threshold operator and model are additional important POCS parameters that need to be determined. Our experiments show no clear preference for either hard or soft thresholding, whereas an exponentially decreasing threshold model is preferred over linear or data-driven ones (e.g., Gao et al., 2010, 2013; Zhang et al., 2015; Wang, 2016; Jicheng et al., 2018). Other adjustable parameters are the maximum and minimum cut-off percentages (⁠$pmax$ and $pmin$⁠) used to calculate the maximum (⁠$τmax$⁠) and minimum threshold limits (⁠$τmin$⁠) based on the maximum sparse coefficient, respectively. For improved noise reduction, the weighting factor $α$ can be reduced accordingly.

When aiming to preserve the highest possible spatial resolution, binning to the “original” trace spacing should be the preferred approach (Figure 14). One drawback of this method is the presence of minor residual artifacts caused by relatively low initial bin coverage and, thus, imperfect interpolation (Figure 14a). Binning the data without sophisticated interpolation using the line spacing (25 m) resulted in inferior quality and a fragmented appearance (Figure 14d). A suitable bin size depends on the average line spacing, which in our case is five times the spatial trace distance along the input 2D profiles. Naturally, closer track lines would allow for binning with higher resolution without sacrificing quality. The optimal bin size for our data set would be approximately 7.5 m × 7.5 m, which proved to be a sufficient compromise of quality and resolution (Figure 14a). When lines with only one predominant orientation are available, it is advisable to increase the bin size perpendicular to these lines accordingly, as our tests show better interpolation results using asymmetrical bins (Figure 14).

We did not attempt to apply time migration to the 2D profiles due to the narrow beam width of the signal emitted by the parametric TOPAS. However, successful migration of nonparametric SBP data has been demonstrated (Baradello, 2014). If applicable to the available data type, 2D migration should be performed immediately prior to 3D binning and subsequent interpolation.

The developed workflow generates comprehensive volumes from sparsely binned 2D profiles (Figure 16) but is not limited to data from hull-mounted parametric SBP systems such as those presented in this study. It also could be applied to data sets acquired using pole-mounted portable systems, towfish instruments, or even autonomous underwater vehicles and uncrewed surface vehicles. The latter options are useful for extremely shallow water settings and hazardous operating conditions. A growing interest in high-resolution systems, particularly for offshore wind farm site surveys, has created an increasing demand for ultra high-resolution (UHR) 3D acquisition systems (Samuel, 2022; PGS, 2023). In addition to dedicated seismic systems such as P-Cable (Planke et al., 2009), recent developments of multitransducer SBPs enable the concurrent acquisition of multiple beams perpendicular to the acquisition heading (Kongsberg Geoacoustics Ltd., 2018; Innomar, 2022). Our workflow could upsample such seismo-acoustic data sets after an initial binning to further increase the spatial resolution. Moreover, we aim to extend the capability of the workflow and perform additional experiments using the full waveform of the TOPAS and UHR seismic data sets.

We present a novel open-source workflow to interpolate sparse 2D subbottom echosounder profiles into a pseudo-3D cube, enabling subsurface features to be interpreted in much greater detail than using only 2D profiles. This workflow is implemented using the Python programming language and its ecosystem of scientific packages. The first workflow stage features the processing of individual 2D profiles to account for vertical offsets by correcting (1) residual static effects, (2) varying tidal elevations, and (3) remaining misties at line intersections. The second stage comprises (1) sparse 3D binning, (2) optional preprocessing steps (gain recovery, amplitude balancing, frequency filtering, resampling, and envelope calculation), and (3) 3D interpolation using the iterative POCS algorithm. Suitable output bin sizes largely depend on the input line spacing and trace distances along the profiles, with bin sizes slightly larger than the average trace distance providing a good compromise between spatial resolution and interpolation quality. The FPOCS algorithm is implemented in the frequency domain operating on frequency slices using the “dask” package for parallel computation, which significantly improves the runtimes of the workflow. In contrast to many other studies, we use the single-scale FFT instead of multi-scale transforms based on curvelets or shearlets. Numerical experiments have shown that the relatively simple FFT is computationally more efficient and generates better interpolation quality for our data set. Currently, our workflow can generate comprehensive 3D volumes from dedicated SBP surveys and profiles acquired concurrently with 3D seismic surveys to provide an additional high-resolution data set. Hence, we are confident that our workflow will be useful for many scientific and industrial legacy data sets. In the future, we aim to further test and extend potential applications such as full-waveform trace interpolation or increasing spatial resolutions of 3D seismic surveys by upsampling through interpolation.

We gratefully acknowledge the Royal Society of New Zealand “Marsden” Fund, grant UOA1022, for funding the research project “Geological champagne: What controls the sudden release of CO₂ at glacial terminations on the Chatham Rise?” This publication is part of the Ph.D. research by F. Warnke that is financially supported by the University of Auckland Doctoral Scholarship. We thank the nautical crew and scientific participants of research voyage TAN2006 with R/V Tangaroa for their support and patience while “mowing the lawn” to acquire the presented TOPAS data set. Ship time was supported by the Ministry of Business, Innovation and Employment (MBIE) and additional funding from GNS Science. Furthermore, we thank the assistant editor, associate editor, and two anonymous reviewers for their constructive comments and feedback that helped to greatly improve the paper.

The open-source workflow is licensed under the GNU General Public License version 3 (GPLv3) and is freely available on GitHub at https://github.com/fwrnke/pseudo-3D-interpolation and additionally archived via Zenodo at https://doi.org/10.5281/zenodo.7647497. Detailed code documentation is available at https://fwrnke.github.io/pseudo-3D-interpolation/. The code and sample data used to generate most of the presented figures are available on GitHub at https://github.com/fwrnke/Warnke-et-al-2023-figures. The TOPAS data are subject to a moratorium of 12 months from the publication date of the paper. Once the embargo expires, the data will be made available in a data repository AND upon request to the corresponding author.

A list of the utilized Python packages for running the workflow is provided in Table B-1.

Biographies and photographs of the authors are not available.