Measurement of sound velocity in oil wells based on fast adaptive median filtering

Measurement of sound velocity in oil wells has long been a challenging industrial issue due to the difficulty in obtaining clear oil pipe coupling waves in strong noise. In this paper, a novel sound velocity measurement method is developed for the dynamic liquid level of oil wells based on fast adaptive median filtering. First, to solve the noise interference problem in the reflected oil pipe coupling wave, a fast adaptive median filtering algorithm is proposed to obtain an accurate oil pipe coupling wave. Then a curve fitting method based on range discrete coefficient is developed to estimate the sound velocity in the tubing-casing annular space of oil wells. The fitting sound velocity can reflect the propagation rule of the real sound velocity. In particular, the sound velocity at any position within the tubing-casing annular space can be accurately calculated by the fitting function. Finally, the proposed method is applied to the dynamic liquid level of an oil well. Experimental results show the effectiveness and favorable accuracy in estimating the sound velocity distribution.


Introduction
Oil is a very important resource for national economies and defense security. The global oil resources are relatively scarce and oil field production capacity is shrinking gradually due to the lack of hydrocarbon resources and low recovery efficiency. Therefore, improving the techniques of oil production and enhancing oil recovery are crucial problems in the long-term strategy for petroleum energy [1][2][3][4][5][6].
There are many factors affecting the recovery of oil wells, including the depth of the dynamic liquid level, the motor speed, and the screw pump. The screw pump oil well must ensure a certain depth of sinking in the process of production in order to have a better effect. Once the dynamic liquid level drops down to the pump's setting depth, the fullness coefficient will be reduced and the friction between the rotor and stator will cause a temperature rise. In this case, the stator will be burned and the well will be damaged. Therefore, it is very important to know the depth of the dynamic liquid level in the oil production process. According to the depth of the dynamic liquid level, engineers can analyze the supply capacity of the well and the relationship between the supply and discharge underground. At the same time, engineers can calculate production parameters such as oil pressure, oil recovery index, and effective permeability by dynamic liquid level. Then they can design the depth of the pump and analyze the factors of energy attenuation. That is very important for designing a reasonable oil extraction technology. Therefore, it is of great significance to accurately measure the dynamic liquid level, which is an effective way to realize oil well stabilization and improve oil recovery efficiency [7][8][9][10].
The echo method, based on the acoustic reflection principle, is a general approach to measure the dynamic liquid level of oil wells. In the echo method, an explosive sound is emitted into the tubing-casing annular space, and the sound is propagated down through the air in the annular space. A schematic of measuring the dynamic liquid level by using the echo method is shown in Figure 1. In the propagation process, the sound pulses are reflected upward through the oil pipe coupling and the dynamic liquid level, and the reflected sound waves are received by the pulse receiver. Then the received sound waves are sent to the filter, amplifier, and signal processing system for computation and processing. The average sound velocity is calculated by the number of reflected oil pipe couplings and the known distance between two oil pipe couplings. Finally, the depth of the dynamic liquid level can be calculated by the average sound velocity and the reflected time. Many researchers have worked on dynamic liquid level measurement, and many different approaches and techniques have been proposed. Meng et al. analyzed the characteristics of the echo curve of the pumping well and designed an echo signal detection circuit, where computers were used to extract the characteristic pulse in the echo signal to calculate the liquid level [11]. The Marr wavelet was employed to filter the interference noise in the echo signal, and accurate oil pipe coupling waves and surface waves can be extracted to measure dynamic liquid levels [12].
Tao et al. adopted wavelet transform to deal with the reflected wave signal [13]. Xiao et al. designed a new dynamic liquid level measuring method to improve the measurement accuracy by upgrading the sound device and the detection method [14]. McCoy et al. proposed a novel method based on sound propagation to measure dynamic liquid level [15]. To prevent the effect of oil foam in high gas oil ratio oil wells, an oil foam formula based on theoretical calculation was proposed by analyzing the flow state in the tubing-casing annular space [16]. Terzic et al. proposed a dynamic liquid level detection method based on ultrasonic sensors and support vector machines (SVMs) [17]. Wang et al. adopted spectral subtraction to deal with background noise in the reflected waves in order to obtain the real reflected signal of the surface and improve the recognition ability for reflected surface waves [18].
It must be noted that the echo method mainly uses the reflected oil pipe coupling wave to estimate the sound velocity and then calculates the depth of the dynamic liquid level using the sound velocity and the time of receiving the reflected surface. Therefore, the sound velocity is directly related to the measurement accuracy of dynamic liquid level. However, sound velocity measurement is very difficult in oil wells due to the nonlinear distribution of the temperature and pressure of air in an oil well. In recent years, many researchers have worked to address this issue. Liu et al. proposed a new method to calculate the sound velocity by indicator diagram [19]. Ren et al. used the infrasonic wave method to measure the dynamic liquid level [20]. In the measurement process, to improve the measurement accuracy, the FFT was used to extract the speed of sound. Liu et al. used the pulse echo method to calculate the sound velocity according to cross-zero detection of the propagation time of the ultrasonic echo [21]. Zhou et al. proposed an average magnitude deference function method to process the coupling's echo wave using many frames, and the time interval of the adjacent coupling's echo was obtained [22]. Then sound velocity was calculated according to the length of the oil pipe.
In the above methods, due to the complicated oil well and the seriously attenuated sound energy, the received signal reflected by the oil pipe coupling and the surface of dynamic liquid level is weak or nonexistent. The received signal often suffers from noise, resulting in oil pipe coupling waves and surface waves not being obvious and the number of visible oil pipe coupling waves being finite. Meanwhile, due to the sound wave having a certain delay in the propagation process, errors are easily produced in the time axis and pulse edge, resulting in errors of sound velocity. Hence, how to accurately estimate the sound velocity in the oil well is the key to measuring the dynamic liquid level, which is directly related to the accuracy of the dynamic liquid level. To address these concerns, we propose a novel sound velocity measurement method for oil wells based on fast adaptive median filtering. First, a fast adaptive median filtering algorithm is developed to deal with the noise in the reflected waves. Furthermore, a curve fitting method based on the range discrete coefficient is presented to estimate the sound velocity in tubing-casing annular space. The fitting formula can reflect the velocity variation of oil wells in a certain area by using the curve fitting method. Hence, accurate sound velocity in oil wells can be obtained. Experimental results show the effectiveness and favorable accuracy in measuring dynamic liquid level.

Fast adaptive median filtering algorithm
Due to the strong noise in the oil well, the acoustic signal often suffers from various kinds of noises in the process of generating and transmitting. Therefore, the echo wave is fuzzy or its characteristics are even submerged. Thus, denoising is very important for the echo wave.

Median filtering algorithm
Median filtering is a simple and very effective filtering method. It can remove and suppress a variety of random noises, so the detailed information of the signal can be well protected. Compared with linear smoothing filters, median filtering can maintain clear waves with the same size of filter window [23][24][25][26][27]. The mid-value reflects the concentration of a set of data, which can report the true value of experimental data better than the mean value. Hence, median filtering has been widely used in the research of miscellaneous subjects.
The experimental steps of traditional median filtering are as follows: • The filter template (a sliding window containing a number of points) is roamed in the signal and the center of the template is coincident with a pixel in the signal.
• Obtain the gray values of each corresponding pixel in the template.
• The gray-scale values are arranged from lowest to highest.
• The intermediate data point of this column is selected and assigned to the pixel of the corresponding template center position.
If there are odd elements in the window, the middle data point is selected as the median value. If there are even elements in the window, the average value of the middle two data points is regarded as the median value. Before processing the signal x(i)(−∞ < i < +∞), we first define a long window with an odd length q = 2m + 1 , where m is a positive integer. Assume that at a certain moment the signal sample in the window is The median filter is a kind of nonlinear signal processing technology based on ranking statistics theory, which can effectively suppress the noise. The basic principle is to replace the value of a point in the digital sequence with the median of the neighborhoods of this point. Hence, the isolated noise can be eliminated. The two-dimensional median filtering output is expressed as where dim[·] represents the median filtering operation, and f (x, y) and g(x, y) are the original waveform and processed waveform, respectively. W is a two-dimensional template, usually using a 3 × 3 , 5 × 5, and 7 × 7 region.
The purpose of denoising is eliminating the noise and protecting the details of the waveform [28]. However, in practical applications, the denoising performance of the traditional median filter is greatly affected by the size of the filter window [29]. There are some contradictions between noise suppression and detail protection. Small filter windows can protect details of the signal better, but the filter ability is limited. With the filter window increased, the noise suppression ability can be enhanced, but it may lead to decrease of detail protection ability and even lost important details. This contradiction is particularly evident for high noise signals.

Fast adaptive median filtering
Considering the existing limitations of the traditional median filtering methods, an improved fast adaptive median filtering method based on the median value is proposed in this section. The basic idea of the proposed algorithm is to determine the noise points and nonnoise points of the signal, and to maintain the gray value for nonnoise points. For the noise point, the size of the filter window is selected adaptively according to the density of the noise, and then fast adaptive median filtering is applied to the noise point. The main steps of the proposed fast adaptive median filtering algorithm are as follows.

Noise detection
The E × F sized echo signal is divided into S subblocks, the z(z = 0, 1, · · · , S − 1) subblock is denoted as B z , the gray value of the detected pixel (b, c) in this subblock is named U (b, c) , and a 5 × 5 sized detection window is formed with this point as the center. The gray value set of all the pixels in the window is  D(b, c) . If the difference between U (b, c) and is determined as the noise point and is marked with Y b,c = 1 . Otherwise, (b, c) is a nonnoise point, marked with Y b,c = 0 . The noise detection method is designed as follows: where H K has a great influence on noise detection and its value is related to the degree of noise interference in the signal. According to experimental results, when the noise is low, we choose a large H K , or else we choose a small H K .

Filter window size
For subblocks B z that have completed the noise detection, the size of the filtering window is determined in real time. The filter window size is determined adaptively by the noise P z in the subblock. P z is the ratio between the number of noise points satisfying Y b,c = 1 and the number of pixels in the subblock. When P z is small, a smaller size filter window should be selected to protect the details of the signal. Otherwise, a larger size window should be selected to enhance the denoising ability. According to this principle, if the square window is selected to filter, the length of window is designed as follows: where w 1 , w 2 , w 3 , and w 4 are four constants satisfying 0 < w 1 < w 2 < w 3 < w 4 < 1 . The following methods are used to determine the four constants' values. First, noises with different frequencies are added to the echo signal. Then we conclude the normalized mean square error (NMSE) and peak signal-to-noise ratio (PSNR) size parameters of the echo signal with different values of four constants. Furthermore, four constants are chosen that satisfy the smaller NMSE and bigger PSNR with different noises. Experimental results show that good denoising results can be obtained with w 1 , w 2 , w 3 , and w 4 as 0.25% , 20% , 40% , and 58% , respectively.

Fast adaptive filtering
After the echo signal is divided into subblocks, the noise detection and the polluted pixel statistics are carried out for each subblock B z . The noise points are filtered out by the fast adaptive filter and the size of the filtering window H z is determined according to the size of the noise. In the proposed fast adaptive filtering algorithm, only a few specific direction pixels' gray values in the window are selected to sort. The window center pixel value is replaced by the weighted medium value of the sorted gray value. Compared to the traditional median filter, where all the sorted pixels' values in the window are used to calculate the window's center pixel value [30], the proposed fast adaptive filtering algorithm can reduce the calculation amount greatly and improve the efficiency. Meanwhile, the fast adaptive median filter has effectiveness in detailed protection of the signal.
After the echo signal is divided into subblocks, the noise detection and the pixel statistics are carried out for each subblock B z , and the size of the filter window H z is determined by the level of noise. Then subblocks are filtered simultaneously so that a great deal of unnecessary waiting time is reduced. Meanwhile, the fast adaptive median filter only processes pixels in a particular direction of echo signal and thus filtering operation time can be significantly reduced and some details of the median filter can be protected.

Sound velocity calculation
In the oil wells, due to the change of oil layer and the nonlinear distribution of the temperature and pressure of air, sound velocity in the tubing-casing annular space cannot be simply considered as a constant [31]. However, it is difficult to obtain information of the oil well at hundreds and thousands of meters underground. Therefore, measuring the sound velocity of the oil well has become a major technical problem, which is directly related to measurement accuracy. In this paper, a curve fitting method based on a large number of measured data in an oil field is presented to obtain the fitting formula that can reflect the sound velocity variation of the oil well. The fitting formula can be used to calculate the sound velocity at any position within the tubing-casing annular space. The average sound velocity in each actual oil well is different, but they are mainly distributed between 310 m/s and 350 m/s. Assume that a single frequency sound pulse is propagated down from the wellhead. Then the propagation time of sound between two adjacent oil pipe couplings is as follows: where T = 1/F s is sampling interval, F s is sampling frequency, and ∆N is the number of sampling intervals.
In theory, the length of each oil pipe coupling is equal to the sound velocity multiplied by the time interval, namely where ν is sound velocity of about 340 m/s and L is 9.6 m, the length of each oil tubing, or namely the distance between two oil pipe couplings.
Assume that the time of emitting sound is t 0 , the time of receiving the first wave crest is t 1 = ∆N 1 /F s , and the time of receiving the n th wave crest is t n = ∆N n /F s , (n = 2, 3, 4, · · · ) . The time interval of the first oil pipe coupling is ∆t 1 = t 1 − t 0 = ∆N 1 /F s . Similarly, the time interval from the (n − 1)th to the n th is ∆t n = t n − t n−1 = (∆N n − ∆N n−1 )/F s , (n = 2, 3, 4, · · · ) . Then the average sound velocity of the first oil pipe coupling of oil tubing is and the average sound velocity of the n th oil pipe coupling is

Curve fitting
The sound velocity in the oil well increases evenly with the rise of pressure. In the actual oil well, the change of the sound velocity with some fluctuation is a nonuniform increase. Therefore, the curve fitting method is used to fit the change of the sound velocity. With the increase of the length of the oil pipe, the sound velocity is fluctuating. However, the tendency of change for sound velocity is sequentially incremented, but the trend is not very obvious. Hence, the fitting function of sound velocity is designed as follows: where α and β are both fit coefficients, ν R is the range discrete coefficient of sound velocity, and ∆f = f s /N is frequency resolution. Sampling frequency is f s = 1000Hz , and sampling points N are 1000, so ∆f = 1 .
Since the increase of sound velocity is slow with the increase of oil pipe couplings, we design the fitting function, where the dependent variable increases with the independent variable, but the increase is smaller.
Meanwhile, we use the range coefficient to reflect the discrete degree of the fitting sound velocity, namely the maximum minus the minimum. The range discrete coefficient indicates the discrete degree of the sound velocity, so the range discrete coefficient ν R is designed as follows: where ν max is the maximum sound velocity, ν min is the minimum sound velocity, and ν mean is the average sound velocity. Since the range discrete coefficient reflects the discrete degree of sound velocity, the sound velocity changes more slowly with a smaller ν R . On the contrary, the sound velocity changes very strongly.

Experiments and discussion
In this paper, the experimental data are obtained from the Dagang Oil Field in China. We collected the sound signal from six different wells from May 1, 2017 to June 1, 2017. The six wells are J1A1, J1B3, J2A2, J3A2, J3B3, and J4A1, and we selected echo sound signals from 29760 datasets, 29760 datasets, 24280 datasets, 19424 datasets, 18784 datasets, and 22630 datasets, respectively.

Fast adaptive median filter experiment
The original reflected oil pipe coupling wave of J3A2 is shown in Figure 2. It is difficult to distinguish the oil pipe coupling wave due to the interference of a large amount of noise. The proposed fast adaptive median filter is used to filter noise. Experimental results are shown in Figure 3. It can be seen that the oil pipe coupling can be clearly distinguished in Figure 3. Experimental results show that the fast adaptive median filtering algorithm has better denoising effects and practical application value.

Sound velocity fitting
According to equation (8), the sound velocity can be obtained by analyzing and calculating the oil pipe coupling's reflected signal. Figure 4 shows the sound velocity curve of six oil wells.   According to equation (10), the range discrete coefficients of the sound velocity are shown in Table 1, where ν max , ν min , and ν mean are the maximum, minimum, and average value of the global data, respectively. The fitted parameter is shown as ν R . It can be seen that the range discrete coefficient is very small, so the sound velocity is changing slowly. According to equation (9), the fitting coefficient of each well can be obtained and the result is shown in Table 2.  Here, α and β are fitting coefficients.
The simulation results show that the fitted parameter ν R · α is much larger than n , so the increasing trend of ν n is slow with the increase of n . The fitting curve of sound velocity is shown in Figure 5. It can be seen that the fitting sound velocity presents an increasing trend, which satisfies the law that the sound velocity will increase with the increase of depth of the oil well.
Due to h = L · n , h is the depth of dynamic liquid level of the oil well. According to equation (9), the relationship between sound velocity and depth of dynamic liquid level can be defined as follows:  In traditional measurement methods of sound velocity, the average of the fitting sound velocity (AFSV) or the coupling measured sound velocity (CMSV) are usually used to calculate the depth of the dynamic liquid level. In this paper, we used the maximum of the fitting sound velocity (MFSV) as the real sound velocity to calculate the depth of the dynamic liquid level. Compared with the AFSV, CMSV [22], and MFSV, the experimental results show that the proposed algorithm offers satisfactory performance.
Further results are shown in Table 3, Table 4, and Table 5, where h r is the practical measured depth, h t = L · n is the theoretical depth, e a is the absolute error, and e r is the relative error. Table 3 demonstrates the superior performance of the fitting sound velocity of our method. Absolute error and relative error both satisfy the requirements of the actual application. Compared with the two calculation methods, the accuracy of using the maximum of the fitting sound velocity as the average sound velocity is much higher than that of the second method. Both the absolute error and relative error are very small, so the results of the measurement values are very close to real values by using the maximum of the fitting sound velocity as the average sound velocity.   We cannot obtain all the oil pipe coupling waves from the wellhead of the oil well to the surface of the dynamic liquid level, due to the seriously attenuated sound energy and the complicated noise in the oil well. Hence, the average of the fitting sound velocity cannot represent the real average sound velocity of the oil, due to the lack of all the reflected oil pipe coupling waves. In this case, the maximum of the fitting sound velocity can represent the real average sound velocity of the oil better, because the real sound velocity increases with the depth of the oil well. Meanwhile, according to equation (8), we can calculate the sound velocity of oil pipe couplings when we obtain more reflected oil pipe coupling waves.

Conclusion
In this paper, we developed a fast adaptive median filtering algorithm to solve the noise problem in reflected waves and a curve fitting method based on range discrete coefficients was employed to estimate the sound velocity in tubing-casing annular space. Experimental results show that the reflected oil pipe coupling signal can be clearly obtained by using fast adaptive median filtering. Then the accurate sound velocity can be extracted easily. The curve fitting function presented in this paper had a good fitting effect on the sound velocity, and the fitting results reflect the law of the real sound velocity in tubing-casing annular space. The proposed method was applied to measure the dynamic liquid level of the oil well, and experimental results demonstrated the effectiveness and favorable accuracy. In future work, we plan to further exploit the relative sound velocity calculation method by improving the interference and adaption of the filter algorithm to enhance effectiveness.