1 Introduction

In recent years, the area of the sea ice in the Arctic has been decreasing with the global warming, and the Arctic shipping lanes are expected to be opened. The Arctic region is rich in mineral resources, accounting for one quarter of the world's oil and natural gas resources. Lots of researches have been conducted in the Arctic region on the oil and gas resources, the hydroacoustic characteristics and compositions of seabed sediments, the distribution and sources of clay minerals, and the climate changes (e.g., Burlin and Shipel'kevich, 2006; Hill and Driscoll, 2008; Zhao et al., 2019; Kim et al., 2020a; Li et al., 2021).

The acquisition of sub-seafloor physical properties and the identification of sub-bottom sediment type and structure are very important for submarine geological hazard surveys. There are mainly two methods to determine the physical properties of the seafloor surface sediments. One method is sampling the sediments by grabbing, box or gravity corer sampling, and then determining the physical properties of the seafloor surface sediments in the laboratory, which is costly and time-consuming (e.g., Meng et al., 2015). The other is the acoustic wave detection method. By recording and interpreting the reflected waves travelling through the sub-seafloor layers, the physical properties of the sub-bottom sediments are obtained. This method is economical and can accurately determine the acoustic parameters, types and physical parameters (e.g., bulk density, porosity, and shear strength) of the seafloor sediments (e.g., Wang and Stewart, 2015; Zheng et al., 2017; Rahili and Cahyono, 2019; Zhen et al., 2023).

Acoustic impedance is a parameter with definite physical significance. It is often used for reservoir characterization and predictive analysis of formation lithology (Wang, 2014), and is the product of rock density and P-wave velocity. Acoustic impedance inversion is simply the transformation of seismic data into pseudoacoustic impedance logs at every trace, with all information in the seismic data retained. Chirp sub-bottom profiler data have been successfully used to quantitatively calculate the physical parameters of shallow seafloor sediments. For example, Schock et al. (1989) calculated the absorption attenuation curve of shallow sedimentary layers using sub-bottom profiler data collected at the Narraganset Bay. LeBlanc et al. (1992) performed acoustic impedance inversion calculations on sub-bottom profiler data from the Narraganset Bay near Rhode Island, USA, and obtained sub-bottom sediment types. Wood and Lindwall (1996) performed full-waveform inversion based on a convolution model on sub-bottom profiler data collected along the coast of Florida to obtain the acoustic impedance of sub-bottom sediment. Bull et al. (1998) calculated the seafloor reflection coefficient based on the amplitude of the primary and secondary reflected waves on the seafloor. Stevenson et al. (2002) calculated the absorption attenuation factor of sedimentary strata from the sub-bottom profiler data using the instantaneous frequency matching method adjusted by a casual attenuation filter. Vardy (2015) extracted the acoustic impedance information of the sedimentary layer from sub-bottom profiler data using a genetic algorithm with the objective function defined by the least square sum of the waveform fitting difference. Liu et al. (2015, 2016, 2020, 2021), Liu and Liu (2015, 2016) and Liu and Zhong (2021) proposed the utilization of developed amplitude variations with offset inversion method to estimate the elastic parameters of the interface between the seawater and seafloor sediment. Chen et al. (2016) used the chirp sub-bottom profiler data to invert seafloor properties in the mud volcano region of the southwestern Dongsha Islands. Kim et al. (2016) explored the buried site of an ancient wooden shipwreck off the west coast of Korea by using chirp sub-bottom profiler data. Li et al. (2017) used the sub-bottom profiler data to analyze the river sediments. The chirp sub-bottom profiler data in the area near the Chukchi Plateau in the Arctic have been widely used for qualitative analysis of the subbottom sediment structures (e.g., Hill et al., 2007; Hill and Driscoll, 2008; Jakobsson et al., 2008; Uchida, 2012; Dove et al., 2014; Kim et al., 2020b). However, there is limited scientific research on quantitatively calculation of the physical properties of sediments. At present, a large amount of sub-bottom profiler data have been collected in the Arctic region (e.g., Zhao et al., 2017; Boggild et al., 2020). Compared to the seabed sampling data, the sub-bottom profile data cover a larger area and are more readily accessible. Hence, the sub-bottom profiler data is conducive to obtaining seabed geotechnical properties more efficiently over a larger area.

In this paper, a sub-bottom sediment acoustic impedance model was established according to the characteristics of the shallow stratigraphic structure at the edge of the deep-sea basin on the west side of the Chukchi Sea, and the effectiveness of the algorithm was verified by iterative inversion. The acoustic impedance inversion (using the time-domain adaptive search matching algorithm) was then performed on the processed chirp sub-bottom profiler data, and the seafloor sediment grain size was determined according to the empirical formula between the acoustic impedance and the sediment grain size.

2 Methods

The time-domain adaptive search matching algorithm first determine the time point m₁ with the maximum peak value in the received waveform record y(n), then takes a waveform signal centered on m₁, which is equal in length to the wavelet, and calculates the cross-correlation value of the waveform signal and the wavelet. When the cross-correlation value is greater than the predefined threshold, this section of signal is considered effective and will be used in the subsequent processes. Otherwise, it is treated as noise and compressed by multipling a very small coefficient. Subsequently, the reflection coefficient at this time point will be calculated by multiplying the peak value M₁ at the time point by a constant scale factor μ, and the reflection coefficient values at other time points are set to 0. Then the remaining waveform record is calculated by subtracting the forward waveform record obtained by the convolution of this reflection coefficient and the wavelet from the received waveform record y(n). The above search steps will be repeated in the remaining waveform record until the iteration stop condition is reached. Through iterations, the reflection coefficient model is updated, the absolute value of the reflected wave amplitude in the remaining waveform records decreases, and the correlation coefficient between the forward waveform records and the initially received waveform records gradually increases.

The specific processes are as follows. Set the maximum value of the sequence y(n) (n = 1, ···, N) as M₁ according to the length L of the incident wavelet x(n) and the waveform shape, identify the reflected wave segment of length L centered on m₁ in the reflected signal sequence y(n), and calculate the correlation coefficient λ of this segment with the wavelet x(n). If the λ value is greater than the given threshold, define:

$ {d_1}(n) = \left\{ {\begin{array}{*{20}{c}} 0&{n \ne {m_1}} \\ {\mu {M_1}}&{n = {m_1}} \end{array}} \right., $

(1)

where μ is a constant scale factor.

Amend the y(n) series to:

$ {y_1}(n) = y(n) - x(n) * {d_1}(n), $

(2)

where '*' denotes the convolution.

Define the d₂(n) sequence as:

$ {d_2}(n) = \left\{ {\begin{array}{*{20}{c}} 0&{n \ne {m_2}} \\ {\mu {M_2}}&{n = {m_2}} \end{array}} \right., $

(3)

where M₂ = y(m₂) is the maximum value of the sequence of y₁(n). Then the current h(n) = d₁(n) + d₂(n), and the p-th sequence of y(n) can be obtained by successive iterations:

$ {y_p}(n) = y(n) - x(n) * \left[ {{d_1}(n) + {d_2}(n) + \cdot \cdot \cdot + {d_p}(n)} \right]. $

(4)

The p-th estimated sequence of h(n) is:

$ h(n) = {d_1}(n) + {d_2}(n) + \cdot \cdot \cdot + {d_p}(n). $

(5)

If the correlation value between the convolved signal h_p(n)*x(n) and the signal y(n) is greater than a given threshold, the p-th estimated sequence h_p(n) is set as the reflection coefficient sequence of the stratum.

In this paper, the value of μ = 0.8 is chosen according to Lin (1997). A μ value greater than 0.8 may make the h_p(n) sequence difficult to converge and hover around the h(n) expectation sequence, while a μ value that is too small makes the convergence slower. Fig. 1 shows the flow chart of the time-domain adaptive search matching algorithm.

When the plane wave is incident normal to the horizontal reflection interface, the reflection coefficient and the acoustic impedance have the following relationship:

$ {Z_{i + 1}} = {Z_i} \times \frac{{1 + {R_i}}}{{1 - {R_i}}}. $

(6)

The recursion formula between the wave impedance value of the (i+1)th layer and the wave impedance value of the first layer is:

$ {Z_{i + 1}} = {Z_0}\prod\nolimits_{k = 1}^{i + 1} {\frac{{1 + {R_k}}}{{1 - {R_k}}}}, $

(7)

where Z₀ is the acoustic impedance of the seafloor surface layer and R_k is the reflection coefficient of the k-th layer. Therefore, once the reflection coefficient of all strata and the acoustic impedance value Z₀ of first layer are accurately known, the acoustic impedance of all strata can be derived from the recursive Eq. (7).

3 Numerical Experiments

According to the characteristics of the shallow sediment structure at the margin of the deep sea basin on the west side of the Chukchi Plateau, a velocity model was established (Fig. 2). From top to bottom, the velocity model consists of 10 layers: sea water layer and sedimentary layers 1 to 9. The density for sea water was set to 1.0 g cm⁻³ and the densities for each sedimentary layer were converted from the velocities using Gardner's empirical formula, i.e., ρ = 0.31·V ^0.25 The velocity, density, and acoustic impedance values are shown in Table 1.

The acoustic impedance model, established based on the velocity model, is shown in Fig. 3. To reduce computation, 37 traces with a trace interval of ca. 28 m were extracted from the ca. 1000-m-long velocity model. The value of the stratigraphic reflection coefficient was calculated from the acoustic impedance model (Fig. 3), and the forward waveform records (Fig. 4) were calculated by convolving the stratigraphic reflection coefficient with the source wavelet (Fig. 5). The source wavelet (Fig. 5) used in the paper is the autocorrelation zero phase Klauder wavelet of a chirp pulse.

The time-domain adaptive search matching algorithm requires input of forward waveform records (Fig. 4) and source wavelet (Fig. 5). And by iteration, the reflection coefficient value is inverted. According to Eq. (7), the acoustic impedance inversion profile is obtained by setting Z₀ equal to the acoustic impedance value of seawater (Fig. 6). The sum of squared errors between the forward waveform records, calculated from the reflection coefficient inversion result and the source wavelet, and the input waveform records is 1.5705×10⁻⁸. In order to visualize the single trace inversion result in the whole acoustic impedance profile, the second synthetic record (Fig. 4) is randomly selected as the input data of the time-domain adaptive search matching algorithm. The output includes reflection coefficient sequence and acoustic impedance values. The reflection coefficients (Fig. 7) are basically consistent with the theoretical model. The acoustic impedance curves obtained from the inversion fit the real acoustic impedance curves very well for each stratum (Fig. 8), with a relative error of less than 2.52×10⁻⁶. Fig. 9 shows that the absolute error between the input reflected waveform data and the inversed reflected waveform record is less than 5.0×10⁻⁵, and the relative error is less than 2.5×10⁻⁶. In this section, we constructed a sedimentary acoustic impedance model of the western deepsea basin margin of the Chukchi Plateau. The Klauder wavelet was selected as the source wavelet, and the theoretical waveform records were obtained through forward calculation. The source wavelet and the theoretical waveform record were used as the input data of the time-domain adaptive search matching algorithm, and the reflection coefficient sequence and acoustic impedance values were obtained through iterations. The sum of squared errors between the theoretical waveform and the waveform calculated from the inversion result is 1.4532×10⁻⁸, indicating the validity of the algorithm.

4 Real Data in the Arctic Ocean

The sub-bottom profiler data for this study were acquired in September 2020 during China's 11th Arctic scientific research cruise aboard R/V Xuelong 2. The experiment was conducted near the Chukchi Plateau area in the Arctic Ocean (Fig. 10). The source excitation parameters are listed in Table 2.

The chirp (Linear FM) seismic wavelet, a signal with instantaneous frequency increasing linearly with time, was used in the sub-bottom profiler measurements on the Chukchi Plateau (Fig. 11). The waveform expression of the seismic wavelet in the time domain is as follows:

$ s(t) = \sin \left({2{\text{π }}\int_0^t {f(t'){\text{d}}t'} } \right) = \sin \left({2{\text{π }}\left({{f_{{\text{start}}}} + \frac{\mu }{2}t} \right)t} \right). $

(8)

The instantaneous frequency f (t) of this signal is linearly related to time and is calculated by the following formula:

$ f(t) = {f_{{\text{start}}}} + \mu \cdot t, $

(9)

where 0 ≤ t ≤ T and T is the time length of the source wavelet. μ is the slope, which is expressed as follows:

$ \mu = \left({{f_{{\text{stop}}}} - {f_{{\text{start}}}}} \right)/T. $

(10)

In this paper, the chirp sub-bottom profiler data were processed using TOPAS software. Firstly, matched filtering is used to improve the resolution. Then, band-pass filtering and adjacent trace stacking are performed to improve the signal-to-noise ratio. The s_seismic2wavelet function provided by the MATLAB toolkit Seislab_10.0301 was used to extract the source wavelet from the data. Subsequently, the processed data and the extracted seismic wavelet were fed into the time-domain adaptive search matching algorithm module to derive the reflection coefficient sequence and acoustic impedance values by iterative inversion. Based on the conversion equation between acoustic impedance and sediment grain size, the seafloor sediment grain size profile was calculated.

5 Results and Discussion 5.1 Acoustic Impedance Profile Result and Evaluation

Eight hundred traces, indicated by the red lines in Fig. 10, were used for the acoustic impedance inversion. Fig. 12 shows the sub-bottom profiler.

After extracting the wavelet, the time domain adaptive search matching algorithm was used to obtain the acoustic impedance profile (Fig. 13). The profile after interpolation in the horizontal direction is shown in Fig. 13, and the lateral tract spacing is 1/8 of the original. Upon comparison of Figs. 12 and 13, the shallow stratigraphic features of the acoustic impedance profile are generally consistent with the stratigraphic pattern on the shallow sub-bottom profile. The average acoustic impedance value of the submarine surface layer in Fig. 13 is 2.5026×10⁶ kg (s m²)⁻¹. The frequency histogram of the acoustic impedance values for the seafloor surface layer is shown in Fig. 14. In the histogram, the values are classified into a total of 80 groups with each group spacing of 0.05×10⁶ kg (s m²)⁻¹. Among all the recorded traces (5161 in total), there are 2640 traces with acoustic impedance values in the range from 2.2×10⁶ to 2.65×10⁶ kg (s m²)⁻¹, accounting for 51.2%. The number of channels with acoustic impedance values in the range from 2.40×10⁶ to 2.45×10⁶ kg (s m²)⁻¹ is the highest with 367 traces. Based on the data from three existing boreholes in the vicinity of this area (S1, S2 and S3 in Fig. 10) (Meng et al., 2015), the wet bulk density values of the seafloor sediments are 1.5, 1.58 and 1.72 g m⁻³, with an average value of 1.64 g m⁻³. The sound speeds in the submarine surface sediment layers are 1425.1, 1441.3 and 1457.1 m s⁻¹, with a mean of 1441.17 m s⁻¹ According to the acoustic impedance formula, the average acoustic impedance value of the submarine surface sediment is 2.3862× 10⁶ kg (s m²)⁻¹, which is well agree with the inversion result obtained in this paper.

To evaluate the inversion result, the distribution of correlation values between the waveform records calculated from each channel inversion result and the input waveform records is shown in Fig. 15 (with a total of 800 traces involved in the inversion before interpolation). Most of the correlation values in Fig. 15 fall within the range between 0.7 and 0.9, with an average of 0.8 marked by the red line. This indicates the strong correlation between the waveform records obtained by inversion and those, and proves that the inversion result is quite reasonable.

5.2 Deep-Sea Basin Grain Size Profiles and Geological Interpretation

According to Meng et al. (2015), the surface sediments in the Chukchi Sea consist mainly of clayey silt. We use the empirical relationship between the Index of Impedance (IOI) of siliceous sediments, the sediment acoustic parameters and physical properties proposed by Richardson and Briggs (2004) (Table 3). The IOI is calculated by multiplying the density of the sediment layer by the ratio of the sound speed in the sediment to the sound speed in the seawater (Richardson and Briggs, 2004). It is independent of the pore water temperature, salinity and pressure, but reflects the intrinsic properties of the medium through which sound waves propagate (Richardson and Briggs, 2004).

The sediment grain size profile (Fig. 16) was calculated according to the empirical relationship provided in Table 3. A comparison with the acoustic impedance profile (Fig. 13) reveals that the magnitude of the acoustic impedance is positively correlated with the sediment grain size.

The mean grain size (θ) of the seafloor surface sediment (Fig. 16) is 7.1498. Fig. 17 shows the frequency histogram of sediment grain size for the submarine strata, classified into 28 groups spaced at intervals of 0.5. The number of samples with grain size values (θ) between 6.5 and 9 is 443, accounting for 68.79% of the total 644 samples, and the group with grain size (θ) between 7.5 and 8 contains the largest number of samples, accounting for 17.7% of the total.

Meng et al. (2015) obtained the grain size values (θ) at three sampling sites (S1, S2, and S3 in Fig. 10) near the study sea area, with values of 7.0, 8.4, and 7.3, respectively (Table 4), which are consistent with our results. According to the mineral grain size classification criteria given in Table 5, the sediment types in the study area are described as very fine silt or clay.

During the China's third Arctic scientific research in 2008, a core, 287.5 cm in length, was collected near the study area (at the site S4 in Fig. 10). A detailed analysis of the samples revealed that the sand content ranged from 0 to 13.73%, with an average value of 1.02% (Chen, 2012). The content of silt is high, ranging from 44.22% to 84.04% with an average value of 60.5%. The clay content ranged from 13.35% to 55.66% with an average value of 38.48%. The sample is silty clay at the top and gradually changes to clay at the bottom, which also shows good agreement with our results.

Base on a cluster analysis of the, grain size and other compositions over more than 4700 core sediment samples from the Chukchi Sea and the adjacent Arctic Ocean, the submarine sediments are divided into 15 categories, with argillaceous silt clay and clay mud being predominant (Kolesnik et al., 2017).

In summary, the seafloor sediment types identified by the method proposed in this paper are largely consistent with the previous sediment classifications in the area, which supports the effectiveness and accuracy of the inversion method.

6 Conclusions

In this paper, an acoustic impedance model was established based on the structural characteristics of the shallow strata at the margin of the deep sea basin on the west side of the Chukchi Plateau, and the forward waveform records were obtained. The acoustic impedance inversion results were obtained after the iterative calculation of the time-domain adaptive search matching algorithm. The misfit between the inversion results and the theoretical model is less than 0.067%. Furthermore, the time domain adaptive search matching algorithm was applied to the inversion of the acoustic impedance over the processed sub-bottom profiler data, which revealed an average acoustic impedance value of 2.5026×10⁶ kg (s m²)⁻¹ and an average grain size value (θ) of 7.1498 for the submarine surface stratum, indicating a predominant distribution of very fine silt.

Acknowledgements

We thank the chief, captain, crew, and science parties of the 11th Chinese National Arctic Research Expedition. This research was supported by the National Key R&D Program of China (No. 2021YFC2801202), the National Natural Science Foundation of China (No. 42076224), and the Fundamental Research Funds for the Central Universities (No. 202262012).