Performance of a modified particle swarm optimization in determining the thickness of the liquid layer outside a well casing pipe

The liquid layer outside a well casing develops from flaws along the cement–casing interface, which have been identified as leakage pathways in wellbores. The thickness of the liquid layer is one of the important parameters for determining cement bond quality. However, it is difficult to accurately determine the thickness of the liquid layer outside a well casing pipe with the traditional ultrasonic pulse echo method. To address this problem, we propose a novel inversion method for determining the thickness of the liquid layer between the casing and the cement sheath. According to test specimens, multi-layered structures are modeled to determine the overall reflected wave at the mud–casing interface. An improved particle swarm optimization algorithm is used to simultaneously inverse the thicknesses of the liquid layer, the casing, and the annular space between the casing and formation. Using synthetic data, this algorithm is more robust and stable when compared with standard particle swarm optimization inversion procedures. In addition, the influences of noise in the reflected wave and deviations of medium parameters on the inversion results of the liquid layer thickness are discussed. Finally, the method was used to determine the thicknesses of the liquid layer in test specimens. Numerical and experimental results demonstrate that the improved particle swarm optimization algorithm provides an effective approach to accurately determine the thickness of the liquid layer outside a well casing pipe.


Introduction
To separate oil-, air-, and water layers, and to prevent leakage and the formation of unwanted interconnecting channels during oil and gas exploration and production, cementing is a commonly used technique in oil wells. More specifically, well cementing is the process of filling cement into the annular space between the wellbore and its casing pipe (Pilkington 1992). The determination of the thickness of the liquid layer outside a well casing pipe is critical to assess the cement bonding performance of different kinds of cements.
The ultrasonic pulse echo method is one of the primary methods for detection of important parameters outside the well casings (Ellis and Singer 2007). The radial measurement technique of this method reduces the effect of the microannulus (Guo et al. 2009) outside the casing pipe and is much less sensitive to tool centralization (Tubman et al. 1984). However, the method has some inherent measurement problems. Numerous studies have tried to eliminate several of these problems. Albert et al. (1988) analyzed the ability of cement bond logs (CBL), ratio bond tool (RBT), and pulse echo tool (PET) logs to detect channels in the casing circumference using a physical model. Pialucha and Cawley (1994) studied the relationship of impedances, layer thickness, and wavelength of the ultrasonic waves employed. Yao et al. (2002Yao et al. ( , 2005 carried out a quantitative inversion of the wave impedance and the thickness of the media outside the casing pipe with a composite reflection coefficient. Che and Qiao (2003) and Ding et al. (2011) analyzed the spectral properties of reflected waves and developed a semi-quantitative detection method for wave impedance and thickness of the casing pipe. Yuning and Xiaoyu (2013) studied the effect of the microannulus on acoustic amplitude. Yuhuan et al. (2014) analyzed the effect of microannulus size on sonic wave propagation. The studies described above revealed the existence of a complex multi-peak-search problem in thickness determination. The above methods overcome this problem by transformations of the objective function. However, the accuracy of the determination methods depends on accurate parameter estimation (e.g., the thickness of the casing pipe). In practice, it is difficult to measure or estimate these parameters accurately.
Particle swarm optimization (PSO) is a populationbased stochastic optimization algorithm (Eberhart et al. 1996;Eberhart and Shi 2001). The algorithm and its different variants have been used across a wide range of applications such as neural network training (Yu et al. 2008), medical diagnosis (Eberhart and Hu 1999), reactive power and voltage control (Yoshida et al. 2000), data clustering (Niknam and Amiri 2010), roundness measurement (Zhu and Chen 2015), machine diagnosis (Guo et al. 2006), data mining (Alatas and Akin 2008), parameter estimation (Donelli et al. 2006), and filter design . In recent years, several limitations of the standard PSO were identified. One of them is a tendency for the swarm to converge prematurely on local optima in complex multi-peak-search problems (Trelea 2003). To solve this problem, a number of researchers proposed variants of the algorithm, which can be classified into several groups: parameter selection, evolution strategy, population topology, and hybrid algorithm.
In this study, we propose a modified particle swarm optimization algorithm to determine the thickness of the liquid layer outside a well casing pipe and discuss the performance of the algorithm after testing it in several numerical and physical models. The objective function of this algorithm is constructed using the spectrum information of the reflected waves at the mud-casing interface. Instead of the thickness of the liquid layer, the thicknesses of the liquid layer, the casing, and the annular space between the casing and formation are inversed simultaneously to eliminate the influences of the parameter estimation errors. In view of the complex multi-peak-search problem in thickness determination and the shortcoming of standard PSO, a modified PSO algorithm is applied to determine the thickness of the liquid layer. Compared to the standard PSO, the static inertia weight is replaced by a dynamic inertia weight, and a dynamic topology and constriction coefficient are used in the modified algorithm.
Physical model and characteristics of the acoustic source Main parameters of the physical model A 2-D schematic diagram of the physical model is shown in Fig. 1. The model consists of five layers with mud, casing, liquid outside the casing, cement sheath, and formation, where d i and Z i (with i = 1,2,3,4,5) denote the layer thicknesses and corresponding wave impedance values of the mud in the casing, the casing, the liquid layer outside the casing, the cement sheath, and the formation separately. The known acoustic parameters of the individual layers are listed in Table 1. Due to the corrosion and deformation of the casing pipes, the thickness throughout the casing pipe is not constant and varies between 7 and 9 mm. The width of the annular space between the casing and the formation d g is the sum of the thickness of the liquid layer outside the casing d 3 and the thickness of the cement sheath d 4 . Because of the eccentricity between the casing and formation, the width of the specific annular space d g varies between 26 and 30 mm. The presence of the rubber layer significantly weakens the reflected waves outside the formation. Thus, they can be neglected in this study.

Characteristics of the acoustic source
At the mud-casing interface, because the acoustic impedance of the casing is much higher than that of the mud, only acoustic waves with certain frequencies can penetrate from the mud into the casing. This specific frequency range is called the casing's resonance and transmission window (Yao et al. 2005). Figure 2 shows the relationship between the reflection coefficient of the mud-casing interface and the frequency of the incident wave. Figure 2 demonstrates that, within the frequency range of the resonance and transmission window (280-420 kHz), the reflection coefficient significantly decreases to a minimum at a frequency of 360 kHz, before increasing again dramatically. Overall, acoustic waves with a frequency of about 360 kHz can penetrate the casing efficiently, which is in accordance with other published literature, for example Che and Qiao (2003) and Yao et al. (2005). Therefore, an acoustic source with an inherent frequency of 360 kHz is selected in our experiments.
Based on the characteristics of the acoustic wave launched from the ultrasonic transmitter, we set up a mathematical model of the acoustic source signal. The function S(t) can be expressed as: In Eq. (1), T s denotes the pulse width of the acoustic signal, with its value set to 16 ls, according to the acoustic signal launched from the transmitter; f 0 denotes the dominant frequency of the signal with a value of 360 kHz. The waveform and frequency spectrum of the acoustic source is shown in Fig. 3.

Numerical simulation and inversion method
Mathematical model of the reflected wave Since the wavelength of the emitted ultrasonic wave is far smaller than the curvature (radius) of the casing pipe, the process of wave propagation through the media in the model can be simplified as the reflection and transmission of a plane wave in a multi-layer medium. Figure 4 shows an N-layered system, of which the left-hand side and righthand side are the two semi-infinite media. The intrinsic acoustic impedance, the wave number of the acoustic wave with a frequency x, and the thickness of each layer are denoted by Z i , k i and, d i , respectively (with i = 1,2,3….N).
Considering both reflection and transmission properties of the plane wave vertically entering a multi-layer planar medium, the expressions of the input impedance of each layer can be written as (Brekhovskikh 2012): where Z ðiÞ r denotes the input impedance of the acoustic wave transmitting from the (i -1)th layer into the ith layer; the intrinsic impedance Z i and the wave number k i of the ith layer can be calculated as Z i ¼ q i c i and k i ¼ x=c i , respectively. In the equations, q i , c i , and x denote the medium density, the Pwave velocity of the acoustic wave, and the frequency of the acoustic wave inside the ith layer, respectively.
The input impedance Z ðnÞ r is equal to the intrinsic impedance Z n of the nth layer since the layer is a semi-infinite medium. The input impedance of the (n -1)th layer can be obtained by Eq. (2). This process is extended to acquire the input impedance of the second layer on the left-hand side of the equation.
The input impedance Z ð2Þ r represents the transformed effect of all impedances of the stacked layers to the right of the first layer.
For the five-layer model shown in Fig. 1, the acoustic wave emitted from the mud layer travels layer by layer and finally is absorbed by the rubber layer. According to Eq. (2), the input impedance of the mud-casing interface can be calculated as: For the acoustic wave, the frequency x, the reflection coefficient RðxÞ, and the reflected wave S r ðxÞ of the mudcasing interface can be calculated as: where Z ð2Þ r is a function of d 2 , d 3 , d 4 , and x according to Eq. (3). As shown in Fig. 1 and Table 1 Combining Eqs. (4) and (3) shows that the reflected wave S r ðxÞ is a function of d 2 , d 3 , d g , and x.

Objective function for the thickness inversion
The reflected acoustic waves within the resonance and transmission window carry information about the property and thickness of the multi-layer medium outside the casing. Therefore, the objective function Jðd 2 ; d 3 ; d g Þ can be constructed from the spectrum of the reflected waves within the window. It is defined as: Þ is the spectrum of the reflected wave acquired from the testing model, and S r ðx; d 2 ; d 3 ; d g Þ is the spectrum of the reflected wave derived from the estimation model. The parameters d Ã 2 , d Ã 3 , and d Ã g are the casing thickness, the thickness of the liquid layer outside the casing, and the annular space width of the testing model, respectively. The angular frequency range of the casing's resonance and transmission window is [x min , x max ]. The three parameters of the testing model are unknown and need to be estimated. The estimation model is a numerical model constructed in accordance with the testing model. The objective function Jðd 2 ; d 3 ; d g Þ reaches a minimum when the parameters d 2 , d 3 , and d g of the estimation model are equal to the corresponding parameters (i.e., d Ã 2 , d Ã 3 , d Ã g ) of the testing model. In this way, we can estimate the values of the casing thickness, the thickness of the liquid layer outside the casing, and the annular Z r (n -1 ) k 1 k 2 k 3 k 4 k n-1 k n d 1 d 2 d 3 d 4 d n-1 d n Fig. 4 Acoustic wave propagation through a N-layered system Z ð2Þ r ¼ space width of the testing model.In the testing model, the thickness of the liquid layer outside the casing, the annular space width, and the casing thickness are 3, 28, and 8 mm, respectively (i.e., d Ã 2 = 8 mm, d Ã 3 = 3 mm, d Ã g = 28 mm). As shown in Fig. 5, the objective function value fluctuates quasi-periodically with the variation of d 3 . Around the global minimum (at d 3 = 3 mm), several local minima can be observed. The effects of these local minima make it difficult to find the minimum of the objective function with linear inverse methods.

Inversion based on an improved particle swarm optimization (PSO) algorithm
For the inversion proposed in this study, the parameters d 2 , d 3 , and d g must be inversed simultaneously. According to the analysis shown in the previous section, there is a complex multi-peak-search problem in the process of inversion. Linear inversion methods (such as the gradient method, the least-square method, etc.) are not effective for this problem unless the initial value of the extremum seeking schemes is selected to be near the global extremum, whereas the standard PSO has the advantage of initial model independence. Some researchers have noted that the standard PSO has a tendency to converge prematurely on local extrema (Trelea 2003). Based on the standard PSO, the static inertia weight is replaced by a dynamic inertia weight, and the dynamic topology and constriction coefficient are used in the modified algorithm. This modified PSO algorithm is called as the KLPSO algorithm.
First, the inertial weight value, which was originally a constant in the PSO algorithm, is now variable in the KLPSO algorithm, i.e., it decreases linearly with the increase of iterations. Second, the global optimum of the particle in the velocity updating equation of the PSO algorithm was replaced with the optimum in the neighborhood of the particle. Third, a constriction coefficient was added to the position of the particle in the updating equation. In the KLPSO algorithm, the velocity, position, and neighborhood of the particles can be updated according to the following equations: where both x and v are three-dimensional vectors. The symbol x denotes the current position of the particle [(x 1 , ; v denotes the current velocity of the particle; k is the iteration number; w is the inertia weight factor; c 1 and c 2 are positive constants and denote the individual and global acceleration coefficients, respectively; b denotes the constriction coefficient; rand 1 and rand 2 are independent uniformly distributed random numbers in the range of [0, 1]; pbest is the individual best position and pnbest is the best position within the neighborhood of a particle; and nRange denotes the number of particles within the neighborhood. The size of the neighborhood increases with the increase of iterations until it includes the entire particle swarm. The initial values of the algorithm parameters are as follows: the number of particles is 20; the acceleration coefficients are c 1 = c 2 = 2.05; the components x 1 , x 2 , x 3 of the position vector are in the range of [7,9], [0, 10], and [26,30], respectively; the constriction coefficient b is 0.729; the inertia factor w is calculated according to Eq. (8), in which CurCnt denotes the current iteration number; the total iteration number LoopCnt is set as 500, and the initial inertia weight w b and the final inertia weight w e are set as 1 and 0.4, respectively.
The KLPSO-algorithm-based inversion procedure is described in more detail below: Step 1 Set the number of particles N = 20 and the maximum number of iterations LoopCnt = 500. Position x and velocity v of each particle are initialized randomly in their respective ranges, and the neighborhood range of the particle is set to nRange = 1. The strategy of uniform distribution was used to initialize the particle position.
Step 2 Evaluate the value of the objective function Jðd 2 ; d 3 ; d g Þ for every particle. This value is then used as the fitness value of this particle.
Step 3 Compare the current fitness value with its own optimum value for every particle. If the current fitness value is superior to its own optimum value, then the current value is used as its individual optimum value and the current position is assigned to the particle's pbest. Define a neighborhood for every particle and determine the local best fitness value and its position for every particle, then assign the position to its pnbest: Step 4 Update the velocity and position of each particle according to Eqs. (6) and (7). Update the inertia weight w according to Eq. (8) and update the neighborhood range of the particles according to Eq. (9).
Step 5 Repeat steps 2 through 5 until a criterion is met (the average of the historical optimum fitness values of all particles is constant or a maximum number of iterations is reached).

Numerical experiments and analysis Stability analysis of the improved PSO algorithm (KLPSO algorithm)
The synthetic reflected wave data set used for the test was derived from a numerical model (Fig. 1), in which parameters d 2 , d 3 , and d g are 8, 3, and 28 mm, respectively. Other model parameters are listed in Table 1. The inversion of the three parameters (i.e., d 2 , d 3 , and d g ) was performed 20 times with the synthetic data set separately by both the KLPSO algorithm and the standard PSO algorithm. Table 2 lists the inversion results of both algorithms. The global extremum is at the point where d 2 = 8 mm, d 3 = 3 mm and d g = 28 mm. As shown in Table 2, the KLPSO algorithm finds the global extremum 19 times, while the standard PSO algorithm finds the global extremum only 6 times. After 20 inversions, the accuracy rate of the KLPSO algorithm is 95%, while that of the standard PSO algorithm is 30%, which shows that the KLPSO algorithm is more stable. Figure 6 shows the evolution of the average optimum fitness value of the objective function for PSO and KLPSO. As the shapes of the curves indicate, the standard PSO algorithm converges quickly and slows down its convergence speed when reaching the optima. The curve for the  Fig. 6 demonstrates that the algorithm has converged prematurely and ended in a local minimum. The curve of the KLPSO shows that the algorithm has a strong ability to avoid the local optima and that it can effectively prevent premature convergence during the evolutionary process. Compared with the standard PSO algorithm, the KLPSO algorithm has better global search ability and significantly enhances the convergence rate and accuracy in this case study.

Effects of the deviations of medium parameters
The deviations of the known medium parameters in the KLPSO algorithm may affect the accuracy of the inversion. As shown in Fig. 7, the effects on the inversion results are induced by the deviations of the acoustic impedance of the mud, the casing pipe, the cement sheath, and the formation when the thickness of the liquid layer outside the casing in the testing model is 3 mm. In Fig. 7, r = (z-z 0 )/z 0 9 100%, which is the relative deviation of the acoustic impedance. The symbol z 0 is the actual value for the acoustic impedance of a certain medium, and z is the value for the impedance of the medium used in the inversion calculation. The symbol d 3 is the inversion result for the thickness of the liquid layer outside the casing. It can be observed that the error of d 3 increases with the increase of the relative deviation of the calculation parameters. The inversion results are greatly affected by the acoustic impedance deviations of the mud and the casing, with maximum relative errors (i.e., the relative error is |d 3 -3|/3 9 100%) of 7 and 6%, respectively.
Both the acoustic impedance deviations of the cement sheath and the formation impose comparatively smaller effects on the inversion results, with maximum relative errors of 0.29 and 0.15%, respectively. Compared with the formation and the cement sheath, the deviations of the acoustic impedance of the mud and the casing have a greater impact on the inversion result.

Effects of noise in reflected waves
The effects of noise on the inversion results using the KLPSO algorithm were explored by adding 10% of white noise to the reflected waves. The inversion values of the three parameters for three different models are tabulated in Table 3. Noise affects the inversion results to a certain degree, but the errors between the inversion values and the actual values are small. From the error analysis, we can conclude that the KLPSO algorithm is able to provide the thickness of the liquid layer outside the casing when the reflected waves contain little noise.   Experiments using a physical model and analysis of the KLPSO algorithm To further test the performance of the proposed algorithm, experiments using a physical model were conducted. Figure 8 shows the experimental device, with the collection system shown on the left and the container holding the physical model shown on the right. Both depth and orientation of the ultrasonic probe in the container can be adjusted using the mechanical device on top of the probe. The physical model in the container has purposely constructed flaws on the casing circumference to produce channels of variable lengths and widths. The model is about 50 cm in height and 60 cm in diameter. The model has two man-made channels embedded in the cement layer. These channels locate at 50°and 230°of the 360°radial surface of the cement-casing interface. Other positions remain good bonding of the cement to the casing. The thicknesses of the liquid layers outside the casing at 16-38 cm height were detected. Figure 9 clearly reveals two man-made channels, and the inversion thickness of the liquid layer does not appear to be much different from the actual thickness. Table 4 lists the inversion results at different sites outside the casing pipe for different thicknesses. When compared to the actual values, the inversion results for the two channels at positions with good cementing quality show small errors. The average errors at for the three different sites are all below 0.18 mm. The average error at the position with good cementing quality may serve as the background error to adjust the inversion results. After the background error is subtracted from the inversion results at the three sites, the adjusted thicknesses are calculated. Their values are very close to the actual values, and the maximum error decreases to 0.027 mm. Hence, we conclude that, using the KLPSO algorithm, the thickness of the channels can be accurately obtained by inversion.

Conclusions
The determination of the thickness of the liquid layer outside a well casing pipe is critical to assess cement bonding performances of different kinds of cements. It is mainly carried out using the traditional ultrasonic pulse echo method, which has difficulty in accurately determining the thickness of the liquid layer.   In this study, a novel inversion method based on a modified PSO was proposed to improve the detection accuracy for the thickness of the liquid layer between the casing and the cement sheath. The method is able to simultaneously inverse the thicknesses of the liquid layer, the casing, and the annular space between the casing and the formation. Compared with traditional approaches, this method does not need to obtain the exact thicknesses of the casing and annular space between casing and formation before inversing the thickness of the liquid layer. Therefore, it can eliminate the influence of parameter errors on the results. Based on the shortcomings of the standard PSO algorithm, the static inertia weight is replaced by a dynamic inertia weight, and the dynamic topology and constriction coefficient are used in the modified algorithm. For the same synthetic data, the accuracy rate of the standard PSO algorithm is 30%, while that of the modified PSO algorithm is 95%, which shows that the modified PSO algorithm is more robust and stable. Using the modified PSO algorithm, we found that the deviations of the acoustic impedance of the mud and the casing have a greater impact on the inversion result than those of the formation and cement sheath. In addition, the new algorithm can achieve good inversion results when the reflected waves contain limited noise. The results of the imaging inversion for a physical model indicate that the modified algorithm can accurately inverse the location and thickness of channels within the cement. This method represents an effective approach to determine the thickness of the liquid layer outside a well casing pipe.