1 Introduction

Ocean bottom nodes (OBNs) have high flexibility, convenient layout and recovery and can obtain a full range of high-fidelity seismic data. These data are conducive to improving the quality of seismic imaging, ensuring the repeatability of 4D exploration, enhancing the effect of reservoir seismic monitoring, and exhibiting great advantages in marine seismic exploration. This approach has gradually become a common observation method for marine seismic exploration (Liu et al., 2017). These OBNs are affected by tides, currents, and other factors when operating for extended periods on the seabed during seismic exploration in shallow sea areas. Throughout the construction period, OBNs may deviate from their initial position and may move several times. Dynamic positioning of OBN is necessary because the conventional secondary positioning technology is not ideal.

At present, the commonly used secondary location methods include arrival wave location and acoustic wave location. Acoustic wave secondary positioning has higher accuracy than the first arrival wave secondary positioning, but it requires additional hardware equipment, is expensive, and is also limited by the recovery rate (Yang et al., 2004; He, 2015). The first arrival wave secondary positioning can not only satisfy the accuracy requirements but is also economical. At present, determination of the secondary location of the first arrival wave mainly involves various methods, such as circle location, near regular tetrahedron method (Yang and Zhang, 2013), refraction first arrival wave location (Zhang, 2005), search method (Ni et al., 2008), variable grid automatic search surface fitting method (Ge, 2015; Di et al., 2016), wave equation method (Huang, 2011), equivalent velocity method (Wang et al., 2020), vector superposition method (Yang et al., 2021), and direct wave and refracted wave joint positioning method (Ma et al., 2023).

At present, all the secondary positioning methods consider the deviation between the actual position of the OBN receivers on the seabed and the dropping position but ignore the multiple movements of the receivers in the construction process. Generally, they only perform secondary positioning, which is inconsistent with the reality of the multiple movements of the actual location of the receivers, resulting in certain errors and affecting the final imaging. In shallow sea OBN exploration, the first arrival seismic wave mainly comprises refracted waves of the seabed and several sets of strata below the seabed. If the first arrival wave includes multilayer refracted waves, then only considering these waves as the refraction of the seabed layer will lead to errors. Therefore, reasonably sorting the shot points corresponding to the refraction layer of the first arrival refraction wave is necessary at different offsets of the receiving point gathers. Subsequently, they are sorted according to different shot times to reasonably assess the relative movements and static time ranges of the receiver points. Multiple-time sorting is carried out using all the shot point information in each period to achieve dynamic positioning. Through this approach, the coordinates of the receiver point corresponding to different shot points are more in line with the actual situation, achieving refined processing and imaging.

In this study, the first arrival time is picked up, and then the sliding horizon of the shot point with different offsets is reasonably divided according to the first-break time and the relative position relationship between the shot point and the receiver. This approach solves the problem regarding the effect of the limited shots of the near offset in shallow sea exploration on the solution accuracy. The corresponding nonlinear equation is established on the basis of the refraction wave positioning principle. The second position of the detection point is identified by calculating the first arrival time. The time segment is optimized according to the time residual, and the positioning is carried out in each period. To overcome the problem of inaccurate first arrival pickup, the use of an optimization algorithm is necessary to reduce the influence of error and determine the optimal solution of the equation (Luo et al., 2022). In this study, the Levenberg-Marquart (L-M) algorithm is preferred to obtain the least squares solution of the secondary location problem through several iterations. Finally, the actual data are used for testing, the L-M algorithm and the moving time of the receivers are used to optimize the final dynamic positioning, and the error between the receiver coordinates and the actual coordinates is small, thereby significantly improving the imaging effect.

2 Refraction Location Method and Optimization of Receiver Moving Time 2.1 Refracted Wave Location Method

In the shallow sea OBN seismic exploration environment, stable refracted waves easily form on the seabed when the seabed topography is relatively flat, the seabed strata velocity exhibits a high-speed layer, and the lateral velocity change is small. The propagation path of the seismic waves follows the minimum Fermat travel time principle. Thus, the direct waves in the first arrival wave are quickly replaced by refracted waves, which are mainly generated and propagated by the high-speed layer on the seabed and below (Xu, 2020). At this time, the coordinates of the receivers can be obtained from the first-break refraction wave received by the receivers (Fig. 1).

In Fig. 1, B indicates the refraction point, and the coordinates of the source and receivers on the sea level are set as S(x_i, y_i, z_i) and R(x, y, z), respectively. The propagation velocity of the seismic wave in water is set as v₁, and the propagation velocity in the seabed is set as v₂. Then, when the seismic wave is excited from point S, the theoretical travel time t_i of the refracted wave generated at point B to the receiver R is obtained as follows:

$ {t}_{i}=\frac{SB}{{v}_{1}}+\frac{AR-AB}{{v}_{2}}=\frac{\sqrt{A{B}^{2}+A{S}^{2}}}{{v}_{1}}+\frac{AR-AB}{{v}_{2}}\text{=}(\frac{\sqrt{A{B}^{2}+A{S}^{2}}}{{v}_{1}}-\frac{AB}{{v}_{2}})+\frac{AR}{{v}_{2}}=\frac{{z}_{i}-z}{{v}_{1}{v}_{2}}\sqrt{{v}_{2}^{2}-{v}_{1}^{2}}+\frac{1}{{v}_{2}}\sqrt{{({x}_{i}-x)}^{2}+{({y}_{i}-y)}^{2}}\text{,} $

(1)

The actual receivers receive the refraction wave, and the corresponding refraction first-break time selected is t_i', X = (x, y, v₁, v₂), i∈[1, N], where N is the number of channels in the common receiver gathers. The error function F(X) is calculated as follows:

$ F(X) = \sum\limits_{i = 1}^N {{{({t_i} - {{t'}_i})}^2}} . $

(2)

F(X) is the sum of residual squares, and the objective is to find x, y, v₁, and v₂, corresponding to the minimum value of F(X) (Gao et al., 2008). Given that the equation is nonlinear, the initial value should be set, and the actual solution should be approximated through iteration. In this study, the nonlinear least squares algorithm L-M is used to solve the secondary location problem.

2.2 Optimization of Receiver Moving Time

To optimize the moving time of the detection point, the L-M algorithm refraction wave positioning method is used for secondary positioning after dividing the subterranean layer, and the coordinates of the detection point after the secondary positioning are obtained.

Using the coordinates of each receiver point, shot point coordinates, seawater velocity, and formation velocity obtained from the preliminary secondary positioning, the theoretical travel time of refracted waves from each shot point to the receiver point is calculated. If the receivers move only once, then accurate coordinates of the receivers can be obtained through the secondary positioning. If the receivers move several times, then the coordinates obtained by the secondary positioning may not be the actual coordinates of the receivers. In this case, deviation is observed between the theoretical first arrival time and the actual pickup time. The difference between the theoretical travel time calculated by the secondary positioning and the actual pickup travel time is determined to obtain the travel time deviation of the refracted wave of each shot. The movement of receivers in the construction process can be judged more accurately and evidently by summing and averaging the travel time of each shot point on the same shot line. The time range and shot range of multiple receiver movements are quickly divided using the deviation, and the optimal movement period of the receivers can be used directly for dynamic positioning.

3 Levenberg-Marquart Iterative Algorithm

The L-M algorithm can overcome the instability caused by the singular and ill-conditioned matrix of the Gauss-Newton method (Xu et al., 2018). At the same time, it is combined with the Gauss-Newton method and the steepest descent method, which is often used to optimize nonlinear least squares problems (Zhu et al., 2016; Wang, 2018).

The nonlinear equations are expressed as follows:

$ \left\{ \begin{array}{l} {f_1}({x_1}, {x_2}, {x_3}, \cdot \cdot \cdot, {x_n}) = 0 \hfill \\ {f_2}({x_1}, {x_2}, {x_3}, \cdot \cdot \cdot, {x_n}) = 0 \hfill \\ ... \hfill \\ {f_m}({x_1}, {x_2}, {x_3}, \cdot \cdot \cdot, {x_n}) = 0 \hfill \\ \end{array} \right. . $

(3)

The vector form is as follows:

$\boldsymbol{F}(X)=0$

$\boldsymbol{X}=\left(x_1, x_2, x_3, \cdots, x_n\right)^{\mathrm{T}}, \boldsymbol{F}(X)=\left(f_1, f_2, f_3, \cdots, f_n\right)^{\mathrm{T}}$.

Eq. (3) can be converted into the following form to solve the nonlinear least squares problem:

$ \mathop {\min }\limits_{x \in {\mathbb{R}^n}} f(x) = \frac{1}{2}\left\| {\boldsymbol{F}(X)} \right\|_2^2 = \frac{1}{2}\sum\limits_{i = 1}^m {{F_i}(X)} . $

(4)

When the minimum value of Eq. (4) is obtained, the gradient and Hesse matrix of the objective function f(x) are obtained as follows:

$ \nabla f(x) = \boldsymbol{J}{(X)^{\text{T}}}\boldsymbol{F}(X), $

(5)

$ {\nabla ^2}f(x) = \boldsymbol{J}{(X)^{\text{T}}}\boldsymbol{J}(X) + \boldsymbol{S}(X), $

(6)

where J(X) is the Jacobian matrix of F(X),

$ \boldsymbol{J}(X) = {(\nabla F{}_1, \nabla {F_2}, \cdot \cdot \cdot, \nabla {F_m})^{\text{T}}} = \left({\begin{array}{*{20}{c}} {\frac{{\partial {F_1}}}{{\partial {x_1}}}}&{...}&{\frac{{\partial {F_1}}}{{\partial {x_n}}}} \\ {...}&{...}&{...} \\ {\frac{{\partial {F_m}}}{{\partial {x_1}}}}&{...}&{\frac{{\partial {F_m}}}{{\partial {x_n}}}} \end{array}} \right), $

(7)

$ \boldsymbol{S}(X) = \sum\limits_{i = 1}^m {{F_i}(X){\nabla ^2}{F_i}(X)} . $

(8)

Some terms of S(X) are omitted for simplicity.

In accordance with the necessary condition for taking an extreme without an optimization problem, when f (x) obtains a minimum at X = (x, y, z)^T, we obtain ∇f = 0. The Gauss-Newton algorithm is used to iterate Eq. (5) to determine the following:

$x^{k+1}=x^k-\left(\boldsymbol{J}_k^{\mathrm{T}} \boldsymbol{J}_k\right)^{-1} \boldsymbol{J}_k^{\mathrm{T}} \boldsymbol{F}^k,$

(9)

where k represents the number of iterations, J_k = J(X^k), and F^k = F(X^k). Iteration step $\boldsymbol{v}^k=-\left(\boldsymbol{J}_k^{\mathrm{T}} \boldsymbol{J}_k\right)^{-1} \boldsymbol{J}_k^{\mathrm{T}} \boldsymbol{F}^k$ is the minimum solution of $ \mathop {\min }\limits_{x \in {\mathbb{R}^n}} \frac{1}{2}\left\| {\boldsymbol{F}({\boldsymbol{X}^k} + {\boldsymbol{J}_k}\boldsymbol{v})} \right\|_2^2 $, which is the

linearized approximation model of the objective function. When the matrix J_k is not a full rank matrix, the approximated formula J_k^TJ_k of the Hesse matrix is a singular matrix, and the Gauss-Newton algorithm is terminated. The L-M algorithm uses regularization term to transform the Hesse matrix into a positive definite matrix J_k^TJ_k + λI. Then, the iteration step is corrected as follows:

$\boldsymbol{s}^k=-\left(\boldsymbol{J}_k^{\mathrm{T}} \boldsymbol{J}_k+\lambda \boldsymbol{I}\right)^{-1} \boldsymbol{J}_k^{\mathrm{T}} \boldsymbol{F}^k ,$

(10)

where I is the identity matrix, parameter λ > 0, and increasing λ can enhance the weight of the identity matrix and improve the number of conditions. Make $ {\lambda }^{k}=\Vert {\boldsymbol{F}}^{k}{\Vert }_{2}^{\delta } $, (1 ≤ δ ≤ 2) so that the value of λ^k changes with the iteration of X value (Fan and Pan, 2006). The updated formula of X is as follows:

$ {\boldsymbol{X}^{k + 1}} = {\boldsymbol{X}^k} + {\boldsymbol{s}^k} . $

(11)

The termination condition of the iterative algorithm is set as follows:

$\left\{\begin{array}{l}\left\|\boldsymbol{g}^k=\boldsymbol{J}_k^T \boldsymbol{F}^k\right\|_2 \leq \varepsilon \\ k>N\end{array}\right.,$

(12)

where the allowable error 0 ≤ ε ≤ 1, and N is the maximum number of iterations.

The specific steps of the secondary positioning method based on the L-M algorithm are as follows:

Step 1. The overdetermined nonlinear equations F(X) = 0 are established according to the shot position data and the first arrival time difference.

Step 2. Parameter δ∈[1, 2], initial point X⁰, tolerance error 0 ≤ ε ≤ 1, and k = 0 are set.

Step 3. g^k=J_k^TF^k is calculated; if ‖g^k‖₂ ≤ ε, then calculation ends. X^k is set as the approximate solution; otherwise,

$ {\lambda }^{k}=\Vert {\boldsymbol{F}}^{k}{\Vert }_{2}^{\delta }, $ $ {\boldsymbol{s}^k} = - {\left({\boldsymbol{J}_k^{\text{T}}{\boldsymbol{J}_k} + {\lambda ^k}\boldsymbol{I}} \right)^{ - 1}}\boldsymbol{J}_k^{\text{T}}{\boldsymbol{F}^k} $ is calculated.

Step 4. X^k+1 = X^k + s^k, k = k + 1 is set.

Step 5. If k > N, then the calculation ends. Otherwise, perform Step 3.

4 Model Test

In this study, a 4000 m × 4000 m × 400 m model is established, in which the water layer measures 0 – 40 m, the first layer below the seabed is 40 – 100 m, the second layer is 100 – 250 m, and the third layer is 250 – 400 m. The receiver is placed at (2030, 1980, 40) (Figs. 2 and 3). The shot interval is 50 m in the x direction over the sea (depth 0), and the shot line spacing is 50 m. The seawater velocity is 1500 m s⁻¹, and the P-wave velocities below the seabed are as follows: 1700 m s⁻¹ in the first layer, 2000 m s⁻¹ in the second layer, and 2500 m s⁻¹ in the third layer. In addition, the receivers are moved several times in the shooting to simulate their multiple movements. The secondary positioning is carried out initially, and then the dynamic positioning trial is performed. The stability of the receivers' movement frequency optimization and the positioning accuracy are verified, considering their differences.

Here, the coordinates of the receivers are moved three times to simulate their drift in the construction process, and then the propagation time is calculated. From the first to the 15th shot lines, the receiver is placed at (2030, 1980, 40). Then, the receiver moves for the first time and is located at (2042, 1988, 40). Subsequently, the 16th to the 38th shot lines are carried out, followed by the 39th to the 55th shot lines at (2065, 2000, 40). Subsequently, the receiver drifts to (2050, 1990, 40), and finally, the 56th to the 81st shot lines are completed.

According to Eqs. (1) and (2), in shallow sea exploration, the propagation time t_i of the seismic wave from the source to the receiver indicates a linear function of the horizontal distance AC between the receiver and the shot point, and the coordinate after secondary positioning satisfies F(X) = 0. The theoretical propagation time of the seismic wave is equal to the actual picking time, that is t_i= t_i'. The mean of offset x_ave and actual propagation time t_ave of all shots within every 50 m spacing is calculated, considering 50 m as spacing and still satisfying Eq. (1). To simplify the calculation, the quotient between the average increment ∆X of the adjacent 50 m offset and the increment ∆t of the adjacent 50 m propagation time is selected. The formation velocity is v = ∆X / ∆t. The propagation layer of seismic waves changes when the v value changes to more than 100 m s⁻¹ (Fig. 4).

The shot point is divided according to the offset. The offsets at 500 and 1200 m are the critical values of seismic wave gliding in the same high-speed layer. The sliding layer of seismic wave is determined; for a detection point, when the location distribution of the shot point and the detection point is uniform and symmetrical, the error resistance becomes robust. Therefore, all the shot point data with an offset of less than 2000 m are selected. The initial position (2000, 2000, 40) is entered. The coordinates of the detection point calculated by the secondary positioning are (2049.02, 1995.55, 41.25), and a certain error with the actual coordinate position is found (Fig. 5). Positioning according to this coordinate may deteriorate the final imaging result; thus, the time should be screened to dynamically locate the receivers.

Through dynamic positioning analysis, the residual time of each shot is calculated, the sum and average of the residual time of all shot points of the same shot line are determined, and the travel time error is obtained (Fig. 6). The absolute magnitude of the change in the average travel time residuals of neighboring gunlines determines the distance that the checkpoints have traveled this value should be judged according to the specific situation and accuracy range. In this study, the value is set to 3 ms according to the actual situation. This approach can determine the movement of the detection point in the adjacent gun line of more than 10 m. If the value exceeds 3 ms, then the receivers have moved. The time is segmented during intervals of relative receiver stability. The receivers are positioned during these stable periods, thereby obtaining results consistent with the actual position (Table 1). This condition is convenient for subsequent research and processing.

5 Application Examples

The field data in this study encompasses 3D marine OBN data. The source is an air gun, which is excited by multiple lines on the sea surface. Its observation system is shown in Fig. 7. The red dots represent shot points, and the blue dots represent the primary location of the receiver points. A total of 64 shot lines with a shot line spacing of approximately 50 m are considered. The sinking depth of the shot point is 6 m. The geophones include the submarine OBNs, and they move several times in one receiving line. The geophone point design distance is 50 m. The sampling rate is 4 ms, and the recording length is 6 s. The water depth in this area is approximately 20 – 80 m, which belongs to the shallow area.

The data from the single shot record (Fig. 8) exhibit a relatively high signal-to-noise ratio, clear reflection wave events, evident refraction waves, and diffraction waves in deep areas.

The actual data were dynamically positioned and calculated to verify the application effect of the method. We obtained the depth by extracting the X and Y coordinates of the shot point and the first positioning coordinates of the corresponding seawater depth from the SPS file or deriving them from the seismic data header to fit the geological survey of the construction area. When the X and Y coordinates of the receivers were obtained accurately, the location point corresponding to the construction area was determined by projection. Combined with the tidal range variation, the depth corresponding to the receiver could be calculated to complete the spatial positioning. According to preliminary calculations, a large number of receivers have moved many times in the X and Y directions but without movements in the Z direction due to the relatively flat sea-bed. Table 2 shows the inversion data of a certain OBN coordinate.

The observation system is redefined after corrections to the coordinates of multiple OBNs. The comparison is made from the arrangement of common receiver gathers, velocity spectrum analysis, and the effect of pre-stack time migration of the first CPM line. Fig. 9 shows the arrangement of common receiver point gathers prior to positioning, after secondary positioning, and after dynamic positioning. Inaccurate coordinates of the receiver points lead to a different arrangement order of the traces in the gather, thereby affecting the continuity of the event. This finding is most evident when viewed from the first arrival time. Prior to positioning, with the increase in offset distance, the first arrival time is chaotic and discontinuous, and the signal-to-noise ratio is poor. After the second positioning, the continuity of the first arrival time is evidently improved. However, in some regional continuities, the result is still not ideal. After dynamic positioning, the first arrival time continuity is the best among the three conditions. Fig. 10 shows a comparison of the velocity spectra. The figure indicates that the part marked in the shallow red circle is dispersed because of the inaccurate position of the receivers, the accuracy of the velocity spectrum is improved after secondary positioning, and the velocity spectrum after dynamic positioning is more convergent with relatively high velocity accuracy.

The profile features in Fig. 11 are compared. The accuracy of the secondary positioning (Fig. 11b) is improved compared with that prior to positioning (Fig. 11a), and the quality of the profile after dynamic positioning (Fig. 11c) is the highest. In particular, the local enlarged view of the red box area can better reflect the advantages of dynamic positioning. After secondary positioning (Fig. 11e), the continuity of the event is better than that prior to positioning (Fig. 11d), but the effect of some areas is unsatisfactory compared with that prior to positioning. This finding indicates that the position of the receiver is still not optimal at this time. After dynamic positioning, the continuity of the event in the migration profile (Fig. 11f) is better and the stratigraphic is more continuous and the horizon is clearer than that prior to and after secondary positioning. The results suggest that the coordinates of the receivers are more accurate after dynamic positioning, and the imaging accuracy is higher.

6 Conclusions

In this study, considering that OBNs in the shallow sea may move several times, based on the refraction positioning principle, the dynamic positioning of receivers is realized by selecting the first arrival of refractive waves, layering the refraction high-speed layer, and conducting time optimization. However, given certain errors in selecting the first arrival and shot point positioning data and because the velocity of the formation is variable, the equations can not be solved directly. The L-M optimization algorithm is used to find the optimal solution for the receivers. Initially, the 3D model data are used for testing, the predetermined position is applied as the initial value input, and the dynamic positioning is carried out through the iterative calculation of the L-M algorithm and the optimization of the time period. These parameters successfully achieve receiver coordinates that are consistent with the actual construction location. Finally, the proposed method is used for dynamic positioning of actual OBN data. The results show that the location of the receivers after dynamic positioning is more accurate than the secondary positioning, and the imaging effect is enhanced.

Acknowledgements

This research is funded by the National Natural Science Foundation of China (No. 42074140), and the Scientific Research and Technology Development Project of China National Petroleum Corporation (No. 2021ZG02). At the same time, we would like to thank the GeoEast Processing System provided by BGP for its support of this study.