SHANG Wenjing ,XUE Wei ,XU Yidong,,MAKAROV Sergey B.,and LI Yingsong

1) College of Information and Communication Engineering,Harbin Engineering University,Harbin 150001, China

2) Key Laboratory of Advanced Marine Communication and Information Technology,Ministry of Industry and Information Technology,Harbin 150001,China

3) Yantai Research Institute and Graduate School of Harbin Engineering University,Yantai 264006,China

4) Institute of Physics,Nanotechnology and Telecommunications,Higher School of Applied Physics and Space Technologies,Peter the Great St.Petersburg Polytechnic University,St.Petersburg 195251,Russia

Abstract The electric inversion technique reconstructs the subsurface medium distribution from acquired data.On the basis of electric inversion,objects buried under the earth or seabed,such as pipelines and unexploded ordnance,are detected and located in a contactless manner.However,the process of accurately reconstructing the shape of the target object is challenging because electric inversion is a nonlinear and ill-posed problem.In this work,we present an inverse multiquadric (IMQ) regularization method based on the level set function for reconstructing buried pipelines.In the case of locating underwater objects,the unknown inversion area is split into two parts,the background and the pipeline with known conductivity.The geometry of the pipeline is represented based on the level set function for achieving a noiseless inversion image.To obtain a binary image,the IMQ is used as the regularization term,which‘pushes’ the level set function away from 0.We also provide an appropriate method to select the bandwidth and regularization parameters for the IMQ regularization term,resulting in reconstructed images with sharp edges.The simulation results and analysis show that the proposed method performs better than classical inversion methods.

Key words inverse problems;level set function;inverse multiquadric regularization method;buried pipeline

1 Introduction

In underwater exploration,artificial facilities such as pipelines and unexploded ordnance may get buried under the seabed because of geological movement or anthropogenic factors.This phenomenon may lead to damaged pipelines or the creation of dangerous zones for vessels (Wu and Guo,2020).Therefore,buried objects need to be detected and located to effectively overcome this problem.Four major noncontact detection techniques,namely,visual-mode techniques,acoustic-mode techniques,magnetic anomaly detection (MAD),and quasistatic electric field recovering method (Adler and Lionheart,2006;Jung and Yun,2014;Simyrdaniset al.,2018;Voet al.,2020),have been developed for detecting buried objects.In general,the cameras used to perform visual detection are mounted on underwater vehicles,such as autonomous underwater vehicles (AUV) and remotely operated vehicles.However,visual detection systems are unable to perform efficiently in turbid environments and are thus unable to detect the target objects buried under the seabed (Chenet al.,2019).Acoustic-based shallow profiler techniques are used to image the target objects buried under the seabed (Tian,2008;Chenet al.,2019).These methods provide a cross-sectional scan of the target area but cannot easily determine the direction of the pipeline intuitively (Tian,2008;Chenet al.,2019).Because most of the pipelines are made of steel,the relative permittivity of the pipeline is larger than seawater and seabed,which may be magnetized by the background magnetic field (Huanget al.,2013,2018;Chenet al.,2019;Jinet al.,2020).Thus,the background magnetic field distibution is distorted by the steel pipeline,resulting in a magnetic anomaly (Huanget al.,2013,2018;Chenet al.,2019).MAD is not an effective technique for detecting nonmagnetic metallic materials.Compared with the visual,acoustic,and MAD technologies,the quasistatic electric field recovering method is not restricted by the water diaphaneity and magnetic properties (Brown,2003).

The quasistatic electric field recovering method is commonly used in locating buried objects,underwater navigation,and geophysics (Wolf-Homeyer,2019).Bazeilleet al.(2017) proposed an electric sense system,which is capable of estimating the pose and geometric properties of nearby objects.Milleret al.(2015) proposed an electric active search trajectory synthesis technique for AUV using nonlinear measurements.Simyrdaniset al.(2018) applied the electric inversion method to recover the 3D image of submerged and buried shipwrecks.

The process of recovering the images of a target based on the electric inversion method poses numerous challenges,such as nonlinearity,ill-conditioning,nonuniqueness,and inverting the matrices with large dimensionality (Sadleir and Fox,2001;Halteret al.,2008;Sunet al.,2019).The electric inversion method also exhibits low resolution and blurred reconstruction in the open domain when appropriate prior information is unavailable (Sunet al.,2019).Various methods have been proposed in the literature to solve the aforementioned issues.The Tikhonov regularization method is the most commonly applied inversion method.This method partially addresses the issues related to ill-conditioning and underdetermining (Sunet al.,2015,2019).The regularization term of the Tikhonov equation depends on the prior information of the spatial characteristics of the target.Thus,the accuracy and stability of the Tikhonov algorithm are strongly affected by the a priori information.Notably,this is one of the key issues in the inversion problems (Borsic and Adler,2012;Ammariet al.,2013).

To achieve a stable inversion and alleviate the image noise,l2-norm regularization terms such as minimum gradient norm and smoothest model (SM) are applied,resulting in a smooth recovered image (Huaet al.,1991;Borsicet al.,2002).Previous studies showed that thel2-norm regularization term-based inversion method outputs images with blurred boundaries.Moreover,spurious oscillations occur during the inversion of a blocky target,indicating that thel2-norm regularization term is unsuitable for recovering targets with clear boundaries,e.g.,pipeline and unexploded ordnance (Portniaguine and Zhdanov,1999;Jafarpouret al.,2009;Ranjanet al.,2018).Notably,the parameters that describe the blocky objects are sparse.Therefore,low-order norm regularization methods are suitable because they tend to reduce redundant information(Ranjanet al.,2018;Shiet al.,2019).Ranjanet al.(2018)proposed a compressed sensing inversion method based on the total variation regularization (l1-norm penalty function).This technique improves the reconstruction resolution as the resulting images exhibit better edge sharpness and reduced image artifacts (Ranjanet al.,2018).The minimum support and minimum gradient support regularization terms are also used in the inversion of blocky objects.These techniques minimize the total domain with the nonzero departure of the model parameters from the given a priori model (Portniaguine and Zhdanov,1999;Huet al.,2019).However,thel0-norm penalty function is non-convex and NP-hard and does not easily converge.Consequently,thel0-norm penalty function is extremely sensitive toward the initial conditions and the settings of the regularization term (Zhdanov and Portniaguine,1999;Huet al.,2019;Utsugi,2019).To ensure the balance between the algorithm’s stability and recovery resolution,hybrid regularization andlp-norm methods are proposed in previous literature.Fournier and Oldenburg (2019) proposed a spatial variable mixedlp-norm based on iteratively reweighted least squares.The inversion result of this technique is controlled by mixedlp-norms comprising multiple terms,each of which is designed to affect the outcome in a specific manner.However,this leads to increased computational complexity.For numerous variables,the settings of the regularization parameters are impractical.Although it shows the full range of solutions,i.e.,0 ≤p≤ 2,and can widen the search,thel1-andl2-norms frequently form the basis in most applications (Fournier and Oldenburg,2019).

Recently,the level set method has emerged as a suitable reconstruction technique for solving the problems related to electrical resistance (ERT) and geology.The level set is a shape-based reconstruction technique (Soleimaniet al.,2006;Kaduet al.,2016,2017;Liuet al.,2017),which concisely describes the material boundaries between the target object and the background,particularly for the anomalies with a constant (known) parameter (Liuet al.,2015;Kaduet al.,2017).Moreover,the level set technique has a strong capability to adapt to the shape changes caused by the merging or splitting of the target (Liuet al.,2015;Kaduet al.,2017).Kaduet al.(2016) proposed a 2D parametric level set reconstruction for reconstructing salt bodies in seismic inversion.The major contribution of this work is the adaptive selection of the narrowband parameter?of the Heaviside function based on the upper bound of the gradient.This setting enables us to achieve a better outcome compared with the original nonlinear least squares problem.

Steel pipelines in the ocean environment can be reasonably approximated as continuous objects with known constant conductivity,which are surrounded by sediment with constant conductivity (Kaduet al.,2016).In this work,we present a method to improve the recovered 3D image of pipelines buried under the seabed.The proposed technique is based on the level set method with anl2-norm inverse multiquadric (IMQ) regularization,and it pushes the estimated parameters to two sides of the level set function.Notably,we expand the level set function using the radial basis function (RBF) to reduce the effective number of parameters.

The remainder of the manuscript is organized as follows.In Section 2,we present a brief introduction of the forward model for quasistatic electric field reconstruction for underwater detection.We also present an overview of the regularization methods specifically used to recover blocky objects.In Section 3,we present the IMQ regularization term-based level set method and solve it using the Gauss-Newton method.In Section 4,we present the numerical results and evaluate the IMQ regularization term-based level set method for different structure buried pipeline targets with high conductivity and low conductivity.We also compare the proposed method with the techniques presented in the literature.Finally,in Section 5,we conclude our work.

2 Theory

2.1 Principles of Electric Field Survey

Electric field inversion is a nonlinear data-fitting technique used to precisely estimate the subsurface properties based on voltage measurements.The electric field inversion process comprises various steps,including 1) the prediction of voltage data by solving the forward model based on an initial conductivity distribution;2) computing the difference between the predicted and measured voltages;and 3) the process of updating the conductivity distribution to improve data fitting.These steps are iteratively repeated until the residual becomes less than a certain threshold(Kaduet al.,2016).To estimate the voltage,a current is injected into the seawater.Consequently,an electric field is generated all over the inversion area.Then,the voltage at different points is measured using the electrode sensors.The electrode sensors have three types of configurations for performing undersea surveys.In the first configuration,all the electrode sensors are located around the surface of the body.This type of configuration is called the closed domain configuration,which is mostly used in ERT for clinical diagnosis (Jung and Yun,2014).In the second configuration,the electrode sensors are located on the surface of the ground.This type of configuration is called the open domain configuration,which is commonly applied in geophysics and undersea target inversion.In the third configuration,all electrode sensors are distributed in the inversion area.This type of configuration is called the body injection configuration,which is commonly employed in geophysics (Constableet al.,1987;Pidliseckyet al.,2007).Compared with the closed and open domain configurations,the body injection configuration usually provides the most accurate recovery.However,the body injection configuration lacks maneuverability.In this work,we use the open domain configuration to recover buried pipes.A detailed survey of the electric field,electrode configurations,and related applications is conducted by Ranjanet al.(2018).

2.2 Forward Model

The conventional quasistatic electric field forward problem in underwater target reconstruction involves the process of solving the potential distribution in region Ω with the boundary condition ?Ω in a given conductivity distribution (Martinset al.,2018).To accomplish the proposed work,we consider a 3D scalar Poisson equation that models the following potential distribution:

Herer∈3R denotes the subsurface position;? denotes the gradient operator;?· denotes the divergence operator;φandσdenote the potential distribution and media conductivity distribution,respectively;andqdenotes the excitation source,which contains the information regarding the position and magnitude of the current injected into the electrodes.

Discrete numerical methods,such as the boundary element method,finite volume method,and finite difference method,are commonly used to deal with the potential distributions in a complex region (Pidliseckyet al.,2007;Shanget al.,2019).A set of linear equations can be derived by discretizing Ω in the continuous function given by Eq.(1) as follows:

where the vectorq,σandφdenote the discretized version ofq,σandφin Eq.(1),respectively;A(σ) denotes the forward model,which is affected by the conductivityσ;andddenotes the prediction data,which are estimated by solving Eq.(2) with the selection matrixQand expressed as follows:

Notably,in the inversion process,the series voltagedis compared with the voltagedobs,which is acquired from the electrodes.Thus,the projection matrixQis used to filter the observed datadfromφ,resulting in the acquired data.

2.3 Inversion Model

The inverse problems are usually nonlinear and ill-posed.Thus,an infinite number of solutions fit the observations equally well.Moreover,the solution may not exist at all.To mitigate these problems,we usually employ regularization techniques.After incorporating the regularization term,the objective function can be expressed as follows:

wherem=lnσdenotes the model parameters,ensuring that the solution is positive (Pidliseckyet al.,2007),andλandRdenote the regularization parameter and regularization operator,respectively.In this work,we assume that the underwater environment can be effectively modeled using a Gaussian channel.The diagonal matrixWdin Eq.(4) represents the weight matrix of the data in each measuring channel for different noise levels (Fournier and Oldenburg,2019;Shanget al.,2020).Thus,Wdfollows a normalized Gaussian distribution,resulting in the data misfit term with a minimum value of 1.Generally,the solution of Eq.(4) is obtained iteratively based on a Newton-like algorithm (Kaduet al.,2016),as follows:

whereγkandHkdenote the step size and an approximation of the Hessian of Φ at thek-th iteration,respectively.The quadratic formulation expressed in Eq.(5) provides a stable solution;however,in several applications,it does not reflect the truth regarding the object that we are trying to recover (Kaduet al.,2017;Xianget al.,2017).For instance,if the parameter is a piecewise constant,thel1-norm-based regularization provides better recovery results.Various techniques have been proposed in the literature to improve the efficiency of algorithms designed to solve nonsmooth problems (Chambolle and Pock,2011;Kaduet al.,2017).In specific scenarios,the object under investigation comprises only two types of materials,which are already known.In these cases,the regularization term is expressed as a nonconvex constraintm∈{m0,m1} (Batenburg and Sijbers,2011;Kaduet al.,2017).Based on this a priori information,we perform the level set reconstruction of the targets,which is more accurate than the traditional reconstruction schemes (Soleimaniet al.,2006).

3 Inverse Problem Using the Level Set Method

3.1 Level Set Method

We assume that there are two areas in the domain Ω with piecewise constant conductivity distribution.The subdomain Ω1denotes the buried pipeline with conductivitym1.The conductivity of the background,i.e.,the seabed,ism0.

Based on the aforementioned assumptions,the inverse problem constitutes finding the edge of Ω1.Thus,the level set functionφis used to formulate a tractable optimization algorithm.The target is located in the area whereφ >0,whereas the background is located in the region whereφ <0.The conductivity distribution is expressed using the level set function as follows:

wherehdenotes the Heaviside function in the interval (0,1) and is mathematically expressed as follows:

However,the ideal Heaviside function is nondifferentiable.Moreover,its gradient cannot be easily computed numerically;therefore,we perform a smooth approximation of the Heaviside function.In recent literature,various functions have been proposed for smooth approximation,including the narrowband,sigmoid,tanh,and arctan functions (Liuet al.,2017).The arctan function is commonly used in various applications.The arctan function is defined as follows:

where the factorμis used to control the gradient ofh.However,to accurately represent the smooth approximation of the Heaviside function,the level set functionφshould tend to ± ∞,which requires a steep gradient ofφnear the image boundary or a largeμ.In Section 3.3,we explain the process of selectingμin detail.To formulate the level set function,we adopt the RBF ?as presented by Aghasiet al.(2011).The RBF is written as follows:

whereαjdenotes the amplitude anddenotes the center of RBF.We rewrite Eq.(10) using the matrix notationφ=Kα,whereα=[α1,α2,···,αN]T,andφ=[φ1,φ2,···,φM]T.The elementkijof the matrixKis defined as follows:

The discretized RBF enables us to represent the unknown conductivity distribution in the form of a vector as follows:

The RBFs are generally used to approximate the smooth multivariate functions,which are commonly applied in the area of high-dimensional data interpolations (Kaduet al.,2016).The most commonly used RBFs are listed in Table 1.

Table 1 Different types of radial basis functions

Notably,the operator (x)+shown in Table 1 returns an output ifx≥ 0;else,it returns 0.The radial behavior of the RBFs is illustrated in Fig.1.The global RBF is highly accurate and often converges exponentially.Therefore,in this work,we use the Gaussian kernel RBF.The inverse problem of recoveringmcorresponds to estimating the distribution ofα.

Fig.1 Variation in the radial basis functions with respect to radius.

3.2 Inverse Multiquadric Regularization Method

The level set method is a binary projection technology,which means that most of the element points are closely populated ‘outside’m0or ‘inside’m1.Furthermore,the level set method aims to prevent the elements from getting trapped in the interval between ‘outside’ and ‘inside’values.Thus,we incorporate the IMQ regularization term based on the RBF.The regularization term acts as a penalty function that severely penalizes the element values near 0.Hence,the objective function expressed in Eq.(4)can be modified as follows:

wheres=Wzφ;Wzdenotes the depth weight,which prevents the anomalous body mainly concentrated on the surface of the inversion area (Pidliseckyet al.,2007);and?denotes the weighting term,which prevents the explosion of the regularization term to infinity when the vector element ofsequals 0.In Fig.2,we illustrate the impact of?on IMQ regularization.For small?,e.g.,?=0.1,the regularization term sets a ‘severe penalty’ within a ‘very narrow’ interval compared with other values.However,the small?may not result in a clear binary image.We conduct further analysis of the impact of?on the numerical simulations presented in this work.

Fig.2 Effect of ? on the inverse multiquadric (IMQ) regularization.

Hence,we use the Gauss-Newton method to minimize the level set objective function expressed in Eq.(13).The gradient can be derived as follows:

where Δd=d(m) -dobsandJdenotes the Jacobian matrix of the forward modeld(m).R1denotes the diagonal matrices expressed as follows:

Meanwhile,we deriveMαusing the chain rule,as follows:

The Hessian matrix is computed based on the gradient of the dual function expressed in Eq.(14),as follows:

whereR2denote the diagonal matrices expressed as follows:

We update the (i+1)-th iteration ofαusing the Gauss-Newton method by combining the Hessian matrixHand the gradient vectorgas follows:

whereηdenotes the step parameter,which is selected based on the algorithm proposed by Pidliseckyet al.(2007).

3.3 Selection of the Bandwidth Parameter

The parameterμin Eq.(9) defines a band in which the Heaviside function is smoothed.In various methods presented in the literature,the selection ofμrelies on the upper and lower bounds of the level set function gradient iteratively (Liuet al.,2017).Contrary to the gradient-based parameter selection method,in this work,we select the value ofμby setting an inexact band of the smoothed Heaviside function.Hence,we can reasonably assume that the inversion process starts from the background parametersφ=φb,h=hband reaches the desired interval whenh≥κorh≤ 1 -κ.Notably,this is equivalent toμ|φ|≥ tan((κ-1/2)π) according to Eq.(9).To provide a clear classification of the background and target pipelines,κgenerally lies in the interval [0.9,1].Thus,we constrain the parameterμaccording toφbas follows:

3.4 Selection of the Regularization Parameter

Generally,the regularization parameterλsignificantly affects the reconstruction results.Several regularization parameter selection methods,such as L-curve,generalized cross-validation,Pareto frontier curve,and discrete Picard condition,are commonly used to solve the inverse problems (Tehraniet al.,2012;Renautet al.,2017;Sunet al.,2019).In this work,the data misfit of the objective function presented in Eq.(13) converges to 1 in an ideal scenario.To obtain meaningful reconstruction results,the model misfit or regularization term should converge to the same order of magnitude with the data misfit.Thus,we propose a simple but effective method to determineλand compare several algorithms,i.e.,the SM,flat model(FM),and modified total variation (mTV) methods.In this work,the level set values of all of the elements areφ=0 ands=0 for the worst-case scenario.The resulting regularization term isλN/2?2,whereNis the number of inversion elements.We restrict the regularization termλN/2?2in the interval [1,10],which results in 2?2/N≤λ≤20?2/N.

In this section,we present the IMQ regularization for recovering the buried pipeline based on the level set method.This method also helps us select the appropriate bandwidth and regularization parameters.Table 2 summarizes the proposed inversion method.In Section 4,we present the performance analysis of the proposed method on the basis of numerical simulations.

Table 2 Process of recovering buried pipelines based on the level set function

4 Simulation Results and Analysis

In this work,we implement the IMQ regularization method using two modified 3D undersea models based on Adam’s synthetic example (Pidliseckyet al.,2007).The first model (model A) is used to reconstruct a curved highconductivity pipeline,and the second (model B) is used to recover branch junctions with low conductivity.We introduce canonical regularization methods,such as the SM,FM,and mTV methods,to compare the proposed method with other techniques presented in the literature.In the case of canonical regularization methods,the regularization parameterλstrongly impacts the reconstruction performance.However,this work aims to determine the effect of the IMQ regularization method used to reconstruct the buried pipeline image.Therefore,the process of selecting the optimalλis beyond the focus of this work.According to work conducted by Borsic and Adler (2012),for a given synthetic model and noise level,we test a series of regularization parameters for determiningλ.This enables us to estimate the reconstructed image that best matches the original image,i.e.,we use the best results provided by each method for comparison (Borsic and Adler,2012).In this work,we perform the simulations using MATLAB.The computer used in this work has an Intel Core i5 3.2 GHz CPU and 8 GB memory.

Model A is illustrated in Fig.3.The figure shows a curved pipeline with high conductivity,i.e.,100 S m-1(the highlight red part).In this simulation,we assume that the background conductivity of the pipeline is 0.02 S m-1(the dark blue part in Fig.3).The dimensions of the model used in this work are 16 m × 16 m × 4.8 m.The model is discretized into a 50 × 50 × 15 mesh grid with a uniform step size of 0.32 m.We consider six emitters and 23 × 23 receivers,which are distributed on the top surface.The blue electrodes are used to acquire the voltage information,and the yellow electrodes serve as source electrodes used to inject current into the seabed.We use the forward model expressed in Eq.(3) to acquire the synthetic data.The acquired data are corrupted with 3% uniformly distributed Gaussian noise.Model B has the same electrode and background conductivity configurations as model A (Fig.4).The conductivity of the buried branch junctions is 10-4S m-1.

Fig.3 Schematic of the numerical computation for reconstructing buried pipeline (Model A),where the background conductivity is 0.02 S m-1.The pipeline is 100 S m-1 (the dark red part).(a),Electrode distribution,where blue electrodes serve as receivers and yellow electrodes serve as emitters injecting the current into the subsurface.(b),x-y cross section of the synthetic model,which is 4 m under the subsurface.

Fig.4 Schematic of the numerical computation for reconstructing the buried branch junction pipeline (Model B),where the background conductivity is 0.02 S m-1.The branch junction is 10-4 S m-1 (the dark blue part).(a),Electrode distribution,where blue electrodes serve as receivers and yellow electrodes serve as emitters injecting the current into the subsurface.(b),x-y cross section of the synthetic model,which is 4 m under the subsurface.

We present an objective and quantitative comparison of various algorithms and parameters.For this purpose,we employ the relative reconstruction error (RRE) previously discussed by Kaduet al.(2016).The RREρis expressed as follows:

wherembdenotes the background parameters distribution with no buried pipeline andmturedenotes the true model parameters distribution.

4.1 Simulation Results

In this work,we use the SM,FM,and mTV methods for comparison.For model A,the reconstruction errors of different methods along with the associated best parameters are listed in Table 3.The initial level set configurations of the IMQ areφb=1 andκ=0.95.According to Eq.(18),we derive the corresponding parameterμ=6.3.The 3D inversion results are shown in Fig.5.Similarly,we present the RRE comparisons in Table 3.

Table 3 Reconstruction errors of different algorithms for model A

The left column in Fig.5 illustrates the inverted 3D entities of the buried pipeline,and the right column in Fig.5 illustrates thex-ycross section of the synthetic model.As shown in Fig.5a,the canonical SM regularization method provides a nearly straight 3D entity image with a small pseudo-abnormal body at the bottom right.This small pseudoabnormal body shows a relatively larger distortion than the true model presented in Fig.3a.Notably,the edge between the background and the pipeline (Fig.5b),where the transition region extends to a large area,is unclear.Moreover,the distribution of the recovered background is uneven,which is not consistent with the actual situation.Fig.5c shows the FM regularization of the 3D image with slight bending properties,introducing a relatively large pseudoabnormal body.Moreover,the background conductivity distribution is not even,as shown in Fig.5d.The 3D recovered image estimated using the mTV regularization method shows the divided structure of the buried pipeline compared with the SM and FM regularization methods,as shown in Fig.5e.We observe the presence of a large number of small abnormal structures in the mTV 3D image,which debase the reconstruction performance.As shown in Fig.5f,the mTV regularization method provides a relatively compact target and even background conductivity distribution.Notably,identifying the edge diffusion of the buried pipeline based on the SM,FM,and mTV regularization methods is difficult.Fig.5g shows a high-quality inversion result.We can easily identify the details and binding properties of the buried pipeline using the IMQ regularization method compared with the SM,FM,and mTV regularization methods.Moreover,the reconstructed background and target conductivity are uniformly distributed with a clear boundary,as shown in Fig.5h.This finding indicates that the IMQ regularization method outperforms the SM,FM,and mTV regularization methods in recovering targets with known conductivity.Table 3 presents a quantitative comparison of the inversion algorithms.Notably,the RRE of the canonical regularization methods,i.e.,SM,FM,and mTV,are greater than 1.3.However,the IMQ regularization based on the level set function provides the best inversion result with the smallest RRE ratio of 0.677.Therefore,this simulation verifies the effectiveness of the proposed IMQ-based inversion algorithm.

Fig.5 (a),(b),image reconstruction results of the buried pipeline using the smoothest model;(c),(d),flat model;(e),(f),modified total variation (mTV);and (g),(h),IMQ regularization methods.(a),(c),(e),and (g),3D entity of the pipeline images;(b),(d),(f),and (h),x-y cross section,which is 4 m under the subsurface.

For model B,the associated best parameters and RRE indices are listed in Table 4.The configurations of the IMQ areφb=1 andκ=0.95,andμ=6.3.The 3D inversion results are shown in Fig.6.We can draw a conclusion similar to Model A from Fig.6 and Table 4,i.e.,the IMQ regularization method exhibits the lowest RRE ratio in recovering buried branch junctions.The SM,FM,and mTV regularization methods provide images with separated structures,and some pseudo-objects are introduced to the inversion results because of which identifying the original image of the branch junctions is difficult,as shown in Figs.6a– 6f.Moreover,the traditional regularization methods could not obtain an evenly distributed background.As shown in Fig.6g,the IMQ regularization method provides images with more details regarding the branch junction pipeline.Fractured parts are detected in the 2D slice of the reconstructed image shown in Fig.6h.However,if by reviewing the 3D image shown in Fig.6g,we can easily observe that the inversion deviation in depth causes this illusion.However,the branch junction structure in model B is more complicated than the curved pipeline structure in model A.Thus,the IMQ regularization method’s RRE ratio of model B is higher than that of model A.In this simulation,the IMQ regularization method based on the level set function also outperforms the SM,FM,and mTV regularization methods in recovering low-conductivity buried objects.

4.2 Effect of the Regularization Term Value

In this work,the regularization termλN/2?2is constrained by the IMQ smoothing parameter?and regularization parameterλ,which significantly affects the resultant reconstructed images.In this subsection,we take model A as an example to analyze the effect of the regularization term value.We fix the initial background parameterφb=1,yielding the parameterμ≈ 6.3 according to Eq.(18).The configurations of?and the correspondingλare presented in Table 5.This results in different values of the regularization termλN/2?2.The resultant reconstructed images of indices 1,2,and 4 in Table 5 are presented in Fig.7.Similarly,the resultant reconstructed images of index 3 in Table 5 are illustrated in Figs.5g and 5h.The IMQ regularization term-based level set method is more capable of recovering the artificial objects buried under the seawater than the canonical regularization methods,i.e.,the SM and FM.All of the 3D images shown in Fig.7 provide a detailed description of the pipeline with clear and uniform background,thus reducing the trailing phenomenon.

Fig.7 Image reconstruction results of the buried pipeline.(a),(b),results obtained with the configuration ? =4 and λ=5 ×10-4;(c),(d),results obtained with the configuration ? =3 and λ=1 × 10-3;(e) and (f),results obtained with the configuration ? =3 and λ=1 × 10-2.

Table 5 Reconstruction errors of different algorithms

As illustrated in Table 5 (indices 2 and 3),the IMQbased regularization method provides inversion results with a low RRE ratio when the value of the regularization termλN/2?2is in the interval [1,10].According to Eq.(4),the data misfit term of the objective function converges to 1,indicating that the model misfit constraints assign more weight at the last stage of inversion.The regularization term of index 1 in Table 5 is smaller than 1,resulting in a higher RRE ratio than indices 2 and 3.Meanwhile,a relatively large regularization term results in a low RRE ratio,as presented in index 4,consequently reducing the adjustment capability of the data misfit term.To achieve efficient and effective reconstruction inversion images,the IMQ regularization term parameters?andλshould follow the restriction1 ≤λN/2?2≤ 10.

5 Conclusions

In this work,we present a reconstruction technique that uses the IMQ regularization term based on the level set method for reconstructing buried pipelines with sharp boundaries.The IMQ regularization term tends to push the distribution parameters away from 0 to a large value,which is equivalent to driving the distribution parametermtoward the binary valuem0orm1.Moreover,the level set scalarμis explicitly derived based on a given background parameterφb.In this work,to obtain reasonable in-version results,the value of the regularization termλN/ 2?2is constrained in the interval [1,10].We use Adam’s synthetic model and modify it by inserting a high-conduc-tivity cured pipeline and a low-conductivity branch junction pipeline to evaluate the proposed method.According to the comparisons presented in this work,the IMQ regularization method provides compact images of the target objects with clearly defined boundaries by reducing the trill phenomenon in thezdirection.In this case,the smallest RRE is 0.6771,thus providing an improved geologic interpretation in terms of edge detection.In our subsequent work,we aim to design methods that can solve undetermined binary distribution inversion problems.

Acknowledgements

This work was supported by the National Natural Science Foundation of China (No.52101383),the Fundamental Research Funds for the Central Universities (No.307 2021CF0802),the Key Laboratory of Advanced Marine Communication and Information Technology,Ministry of Industry and Information Technology (No.AMCIT2101-02),the Sino-Russian Cooperation Fund of Harbin Engineering University (No.2021HEUCRF006),the Ministry of Science and Higher Education of the Russian Federation(No.075-15-2020-934),and the International Science &Technology Cooperation Program of China (No.2014DF R10240).