Research on spectral method of lost circulation layer located by transient pressure wave

Leakage in the wellbore annulus during drilling operations can affect normal production operations, resulting in a severe waste of resources and economic loss, so it is crucial to adopt a fast and effective leak identification method for subsequent plugging operations. For the complex problem of judging the location of the leakage layer, we proposed the method of excitation pressure wave to identify the location of the leakage layer. By analyzing the transfer of pressure waves within the annular pipe system and the pressure head response spectrum, the leak's location is identified based on the location of the resonance point and the change in resonance amplitude. The pressure wave signal contains too much noise. The variational mode decomposition (VMD) algorithm and the Hilbert joint spectrum were used to extract the main frequency components to reconstruct the signal to achieve the denoising effect. On this basis, the reconstructed signal is processed by fast Fourier transform (FFT) to obtain the pressure wave response spectrum, analyze the frequency domain features, and then determine the location of the leakage layer. The experimental results verify that: ① When a leak occurs in the wellbore annulus, the pressure wave will generate additional resonance points in the frequency domain due to the presence of the leak point. ② The combination of the VMD algorithm, FFT, and Hilbert joint spectrum can effectively remove the noise of the pressure wave signal. ③ The method effectively avoids the difficulty of identifying negative pressure waves in the time domain analysis of pressure wave signals. It can effectively locate the leaky layer in the frequency domain analysis. It is concluded that the principle of the method is feasible and has practical significance for field application.


| INTRODUCTION
When drilling deep formations, leakage occurs due to frequently drilling porous and highly fractured formations. [1][2][3] Once the wellbore leaks, the loss is also huge. For years, leakage prevention and plugging have been the research focus in drilling engineering. [4][5][6][7][8] The key to successful plugging is to quickly and accurately determine the location of the lost circulation layer, which is conducive to shortening the plugging time and reducing the total cost of drilling.
There are three broad categories of methods to determine the location of the leaky layer: direct observation, machine learning prediction, and instrumental testing. The direct observation method is to observe the response in the drilling construction process and to analyze all the logging and logging information to determine the location of the leaky layer, which has uncertainty and low reliability. Machine learning models can predict lost circulation, such as lost circulation type and estimated lost circulation amount, and can provide suggestions for preventing lost circulation and making remedial decisions. 9 The most common machine learning methods for prediction are artificial neural networks (ANNs) 10 and support vector machines (SVMs). 11 Abbas et al. 12 built a model based on drilling data using ANNs and SVMs that can quickly and efficiently predict well leaks. Sabah has developed an intelligent prediction model including decision tree, 13 adaptive neuro-fuzzy inference systems, 14 ANNs, and hybrid ANN algorithms, 15 which can quantitatively predict lost circulation. 7 In practical application, these methods need a large amount of data to train a good machine learning model for establishing a lost circulation prediction model. The machine learning prediction method can only provide predictions and suggestions for possible leakage. It cannot directly determine the leakage horizon, so it is challenging to promote its use. Compared with the direct observation and machine learning prediction methods, the instrumental measurement method is more reliable and effective through instrumental measurement.
The instrument testing method is to use logging instruments, such as electronic flowmeter, temperature measuring instrument, pressure measuring instrument, radioactive tracer, resistance temperature detector (RTD) measuring method, and other professional instruments, using sensors to detect the flow rate, pressure, temperature, and other parameters in the well, and then to collect and process the data through the circuit. According to the characteristics of the data variation in detecting different subsurface locations, the actual measurement data is used to locate the location of the leaky layer by integrating the specific construction situation. Radiotracer measurements, 16 thermal resistance measurements, pressure models, [17][18][19][20] and well temperature measurements 21 are widely used in instrumental testing methods. The radioactive tracer measurement method involves pumping drilling fluid containing radioactive material into the well and determining the location of the leak depending on the difference in the natural gamma curve of the radioactive tracer. The RTD measurement method is to put the RTD instrument into the predicted leakage layer downhole, record the resistance value at this time, and then pump the drilling fluid into the well. If the resistance value changes after pumping, the leakage layer is below the RTD instrument. If the resistance value remains the same, then the leakage layer is above the instrument, this method has a certain blindness, and the leakage measurement is inefficient. In 2016, Bogdan et al. 22 proposed high-frequency pressure monitoring (HFPM) for fracturing detection. HFPM is a new noninvasive method implemented by advanced signal processing and interpretation algorithms for monitoring multistage operations and making real-time decisions, identifying downhole conditions and leak locations based on pressure variation profiles collected at the wellhead. In 2016, Chen and colleagues [23][24][25] established the first fluid heat conduction model under the condition of leakage circulation (Chen model for short), as well as proposed to use of microchips in the drilling fluid to extract good data while maintaining mud circulation, predicting the location of loss when a leak occurs based on the established wellbore temperature model. However, the Chen model ignored the effect of the heat source term on the wellbore temperature distribution in the build-up process. In 2019, Wu established a temperature field model of a straight wellbore under leakage cycle conditions, which can determine the location of the leakage layer based on the distribution curve of temperature gradient in the annulus. Compared with the Chen model, the temperature calculation results of the proposed model are closer to the measured temperature. 26 However, the method requires accurate measurement of the wellbore temperature distribution under leakage conditions and the development of supporting wellbore fluid temperature measurement instruments, so there are limitations in its practical application.
The annulus system has the characteristics of uniform, vertical, fluid-like pipelines. Therefore, the problem of lost circulation identification can be studied in combination with pipeline leakage detection methods. In 2001, Mpesha et al. 27 first proposed the impulse response detection method and the transfer matrix method for pipeline leak detection. A steady oscillating water hammer wave is generated by periodically opening and closing the valve. Considering the stable oscillating flow, a frequency response analysis is conducted on the system with leakage to identify the leakage location. This method only requires adding valves and pressure sensors at the equipment outlet under test, characterized by simple construction and fast and effective response. Many scholars have continued and studied the method from the time domain, [28][29][30] the frequency domain, [31][32][33][34][35] and signal processing. 36,37 In 2003, Ferrante applied the impulse response method to solve the control equations for transient flow in pressure pipelines directly in the frequency domain and through the expressions for the pressure spectrum at the downstream end of a single pipe system to analyze the pipeline leakage system. 38 This method shows that the pressure head spectrum is the product of two terms: the former is related only to the maneuvering characteristics, and the latter depends only on the system characteristics. From a technical point of view and considering the field operation, the method of transient excitation proposed by Ferrante et al. is more maneuverable and reliable than the method of steadystate oscillation proposed by Mpesha et al. Lee et al. 39 analyzed the transient trajectories generated by fluid transients injected into the pipeline in the frequency domain and proposed a method for pipeline leak detection based on the frequency response of the pipeline, and determined the leakage area by comparing the relative magnitude between the peaks in the frequency response graph. Lee extracted the system frequency response function using the valve closure method and compared the corresponding time domain effects for other working conditions such as friction and blockage. 40 Ghazali et al. 35 analyzed the pressure transient leakage signal of the water supply network and analyzed the transient frequency characteristics of the signal by signal processing using the Hilbert transform, cepstrum, and other methods. Scola 41 proposed a fault identification method to detect blocked and leaking pipes and validated it with an example. Capponi et al. 42 analyzed the pipeline time and frequency domain model and proposed using correction factors to improve the model effect. To sum up, the pulse response method does not require complex measuring instruments in pipeline leakage detection, and it can detect and verify the leakage area from time and frequency aspects, which can narrow the detection error range, so it has the advantage of being fast, effective, and accurate. In this paper, we propose a new method to locate the loss zone location during drilling by the impulse response of the transient pressure wave while maintaining drilling fluid circulation.

| METHOD AND MATHEMATICAL MODEL
Based on the induced response, the impulse response method and the convolution theorem are applied to derive the expression of the pressure response in the frequency domain at the outlet of the hollow loop end. From the expression, it is known that the transfer function is only related to the system's characteristics, and the Fourier transform of the flow variable is only related to the operating characteristics. The pressure frequency domain response and the Fourier transform of the flow variable are linear, according to which the critical information of the pipeline leakage point can be obtained in theory. The presence of a leak will cause a change in the resonance characteristics (system transfer function) of the wellbore annular system. Under the premise of constant operating characteristics, when pressure fluctuations are excited at the end of the wellbore annular, additional resonance points will be generated in the frequency domain due to the presence of the leak. The leak's location and the leak flow rate will significantly affect the location of the resonance points and the magnitude of the resonance amplitude. The location of the resonance points and the magnitude of the resonance amplitude between the measured pressure head spectrum and the theoretical response spectrum are analyzed to determine the leak's location in the wellbore annulus.

| Transient pressure and flow transfer model in the annulus
A vertical wellbore model of a leaking wellbore annulus is shown in Figure 1 to simulate the flow and pressure transfer process of drilling fluid when a leak occurs in the wellbore annulus. When a leak exists in the wellbore annulus, the leak divides the wellbore annulus system into two parts, annulus 1 and annulus 2, with the length of annulus 1 being L L L 1 = − 2, from the bottom of the annulus to the leak; and the length of annulus 2 being L2, from the leak to the valve. A valve is added at the choke manifold to generate an excitation pressure wave. When the drilling fluid is circulated back out, the valve is partially closed instantaneously and quickly, and the excitation pressure wave is generated due to the water strike effect.
The transfer matrix method is very effective for the temporary pipeline flow problem caused by the excitation wave. According to the distributed parameter theory, the terminal pressure and flow of the annulus with length can be expressed by the initial pressure and flow.
The field matrix equation of a one-dimensional pipe can be expressed as follows 27 : where H D represents the terminal pressure; Q D represents the terminal flow; H U represents the initial pressure; Q U represents the initial flow; γ Cs Ls R = ( + ) is the propagation constant; Z γ Cs = / C represents characteristic impedance; L gA = 1/ represents impedance per unit length; s iω = represents Laplace variable; C gA a = / 2 represents capacitive reactance per unit length; R fq gDA = /(2 ) 2 2 represents impedance per unit length; g represents gravitational acceleration; A represents the cross-sectional area of the annulus; a represents pressure wave velocity; D represents the inner diameter of the annulus; f represents Darcy-Weisbach friction coefficient.
. The relationship equation holds, Z L represents the dynamic impedance at the loss zone. Therefore, the transfer matrix equation at the loss zone is shown in Equation (2) as The total transfer matrix equation of the wellbore annulus can be obtained by combining Equations (1) and (2) as shown in Equation (3) as The impedance at any point within the wellbore annulus is defined as follows: When studying the leakage of the wellbore annulus system, as shown in Figure 1, since the beginning of the annulus is a constant pressure source, the fluctuating pressure at the inlet is H = 0 U1 , and the impedance at the inlet of the annulus is . The cross-sectional area of the annulus remains the same. Therefore, γ γ γ = = 1 2 and (3), the output impedance of the annular is shown in Equation (5) when there is no leakage in the annulus; when there is leakage in the annulus, the output impedance of the annulus is shown in Equation (6). Equations (5) and (6) are defined as the transfer functions of the wellbore annulus system 2.2 | Spectral analysis of transfer function in the annulus

| Spectrum analysis of transfer function for the annulus without the loss zone
When friction exists in the wellbore annulus, the propagation constant is expanded in power series as follows: Neglecting the higher order terms of R and taking the first two terms. The transfer function of the annulus is as follows: where represents the dimensionless angular frequency.

| Spectrum analysis of transfer function in the annulus with the loss zone
When there is a loss zone in the wellbore annulus, the annulus is divided into two parts for analysis. According to the overall transfer function of the annulus to analyze the characteristics of the annulus system, ignoring the pipe friction, γ γ hold, then the output impedance of the annulus is as follows: According to Equation (9), it can be concluded that when the wellbore annulus leakage is minimal, Z L tends to infinity, then Z Z C L tends to zero; when the wellbore annulus leakage is substantial, Z L tends to zero, then Z Z C L tends to infinity.

| Frequency response of flow variable
The beginning of the annulus is a constant pressure source. When the outlet of the annulus is imposed on the flow variable due to the momentary partial closure of the valve, the relative flow will cause a similar step signal change as follows: When the flow variable is given in the time domain or accurately controlled during valve operation, its Fourier transform can be obtained to obtain its frequency response diagram. The Fourier transform of the variable flow function is as follows: The frequency domain amplitude mode function is as follows: where q Δ ′ represents the steady-state flow rate of the annulus minus the average flow rate after the valve is closed; T = C l a 2 represents the annulus half-cycle. The pressure at the outlet of the annulus is constant. When the pulse flow is imposed on the outlet of the annulus due to the immediate partial closure of the valve, the Fourier transform of the unit pulse function is as follows: According to Equation (14), the unit pulse is the sum of infinite sinusoidal excitations. Taking the exit of the wellbore annulus as the starting point, find the impulse pressure response and impulse flow response at any position in the wellbore annulus: When x = 0, the impulse pressure response and the impulse flow response at the valve are as follows: Equation (15) is called the impulse response function. Set the flow variable imposed at the valve, and the flow variable function of the instantaneous partial closure of the valve is the Heaviside function. The Fourier transform convolution theorem can be used to calculate the pressure head function relative to the impulse pressure response function and the flow variable: Equation (16) is the analysis of the pressure head function in the time domain. Apply the nature of the convolution theorem, the Fourier transform of the pressure head function can be obtained as the product of the Fourier transform of the impulse pressure response, the Fourier transform of the flow variable, and the pressure head in the frequency domain is derived as follows: Equation (17) shows that the pressure head spectrum consists of two components for a valve operation with an immediate partial closure: the Fourier transform of the system transfer function and the Fourier transform of the flow variable. The former is only related to the system's characteristics, the latter is only associated with the valve operating characteristics, and the pressure head spectrum and the Fourier transform of the flow variable are linear.

| SIMULATION ANALYSIS
The wellbore annulus system is simulated according to the transfer matrix equation and the transfer function of the annulus, and the parameters in the simulation process are shown in Table 1. When the annulus is without the loss zone, the transfer function spectrum of the annulus system is obtained from Equations (1) and (4), as shown in Figure 2.
In case of leakage of the annulus, the loss zone at 1/2 and 2/3 of the wellbore annulus with a total length of 1000 m shall be simulated with the wellhead as the starting point according to the location of the loss zone of the annulus. The spectrum of the system transfer function when the annulus leaks are obtained from the leak point transfer matrix (2) and the total transfer matrix (3), as shown in Figure 3.
The transfer function is only related to the system characteristics, and the existence of the loss zone will significantly affect the system characteristics. It can be seen from Figures 2 and 3 that the frequency spectrum of Friction coefficient f 0.5 the system transfer function will be substantially affected when the wellbore annulus leaks. The impulse flow is imposed at the outlet of the annulus due to the rapid closing of the valve, causing changes similar to step signals. The impulse flow is used as the input to the annulus system, and the impulse pressure response and flow response are used as the output of the annulus system. The spectrum of flow variables is obtained according to Equation (11), as shown in Figure 4. According to Equation (17), the Fourier transform of the pressure head spectrum and the flow variable can be derived as a linear relationship, and the pressure head spectrum is shown in Figure 5 when leakage occurs at 1/ 2 and 2/3 of the distance from the wellhead.
In Figure 5, comparing the spectrum of the pressure head function when the leakage is located at 1/2 and 2/3 of the wellbore annulus, respectively, starting at the wellhead, it can be obtained that the different locations of the leakage significantly change the spectrum of the pressure head function with the same operating characteristics of the control valve closure. Therefore, for leaks at various sites in the wellbore annulus, the resonance points and resonance amplitudes of the pressure waves in the frequency domain are analyzed to locate the location of leaks in the wellbore annulus. Figure 6 shows the wellbore annulus test system containing the central control computer, water tank, drilling pump, leakage valve, choke manifold valve and pressure sensor,  and high-frequency acquisition card. The experimental wellbore was built as a straight well with the parameters shown in Table 2. The two leakage valves are located at 18 and 24 m of the straight wellbore with a total length of 36 m. A valve is set at the choke manifold at the top of the wellbore to generate the transient wave. The high-frequency acquisition card can store the data collected by the pressure sensor and use the collected data to analyze the change of the pressure wave before and after the evolution of the annulus characteristics. The test steps are as follows: inject water into the water tank, open the drilling pump for circulation, and keep the valve at the choke manifold fully open. After 8 min of average circulation, open the leakage valve on the wellbore to create leakage conditions. After the fluid in the wellbore is stabilized again, record the steady leakage amount. The central control computer controls the valve at the throttle manifold to close momentarily to a specific opening degree and then open quickly. Highfrequency acquisition cards capture and record pressure sensor data from the start of excitation until the pressure waveform disappears.

| Experimental methods
When leakage occurred at 1/2 and 2/3 of the wellbore annulus, the time domain signal of the pressure waveform collected by the HF acquisition card is shown in Figure 7, from which it can be obtained that the pressure waveform jumped significantly after the valve was partially closed quickly. To further get information on the annulus leakage location from the frequency spectrum of the pressure wave signal, it is necessary to convert the pressure wave signal from the time domain signal to the frequency domain signal. The fast Fourier transform (FFT) is a transform method for analyzing signals in the time-frequency domain. It has an extensive range of applications in signal analysis, with simple and intuitive features. In 1965, Cooley and Tukey 43 proposed a fast algorithm for discrete Fourier transform (DFT), namely, FFT, which significantly reduced the computational complexity of the original DFT solution. The spectrum of the pressure wave signal after an FFT is shown in Figure 8. Both the pressure wave spectrum in Figure 8 and the pressure waveform in Figure 7 shows that the acquired signal contains too much noise. Due to the interference caused by the pump and other equipment to the pressure sensor during the experiment, the acquired signal contains too much noise, making extracting helpful information directly challenging. Hence, the received signal needs to be denoised.

| VMD algorithm for signal denoising
For signal denoising, the wavelet transform method, empirical mode decomposition (EMD), sparse decomposition, variational mode decomposition (VMD), and the improved algorithm based on them are the denoising methods adopted by many researchers at present. [44][45][46][47][48] In all of the above methods, there is a certain degree of denoising ability for the signal, but there are also shortcomings. The wavelet transform method has strong applicability to nonsmooth nonlinear signals. Still, the denoising effect is affected by the characteristics of the signal itself, and it is also essential for the selection of the threshold and the quantization of its coefficients. Empirical mode decomposition (EMD) is a timefrequency domain analysis method that removes the Fourier transform's limitations and can be decomposed  adaptively according to the characteristics of the signal itself, but it lacks strict mathematical theory support. EMD also has the phenomenon of mode aliasing when the characteristic signal scale changes intermittently. 49,50 Most sparse decomposition methods are aimed at a single signal, and the required atomic library must be too complete. The VMD algorithm is an adaptive and fully non-recursive signal denoising method proposed by Dragomiretskiy and Zossoin 2014, 51 which mainly constructs and solves the variational problem and then transforms the variational problem into a nonvariational problem to solve, which can decompose a complex nonsmooth nonlinear signal into multiple band-limited intrinsic mode function (BIMF) and a residual signal. The VMD algorithm effectively reduces the problems of modal aliasing, boundary effects, and overenveloping. [50][51][52] It is more sensitive to low-frequency features and more suitable for processing pressure wave signals whose helpful information is contained in the lowfrequency band, 53,54 so it is more appropriate to use this method to denoise the acquired pressure wave signals.

| Decomposition principle of the VMD algorithm
VMD algorithm decomposes the original signal into multiple BIMF and the sum of residual. The number of decompositions and the second penalty operator are critical parameters of the VMD algorithm, which must be set in advance. The IMF is redefined in the VMD algorithm as an amplitude modulation-frequency modulation signal with finite bandwidth, and the following equation defines the decomposed BIMF: The K modal function obtained from the original signal decomposition has its central frequencies, based on which the constrained variational model is solved to determine the bandwidth of each BIMF, and the corresponding variational constrained model is as follows: The constrained variational problem is transformed into an unconstrained variational problem by introducing the quadratic penalty operator and the Lagrange multiplicative operator to solve the constrained model and obtain the augmented Lagrange expressions as follows: The alternate direction method of the multiplicative operator is used to iteratively update each component and its frequency center to obtain the optimal solution to F I G U R E 9 Time-frequency diagram of each band-limited intrinsic mode function (BIMF) after pressure wave decomposition when the leakage position is at 1/2 of the annulus. the initial problem. The expressions for the modal functions and the center frequencies are as follows: For example, when solving the model in Equation (23), the number of decomposition BIMFs and convergence tolerance can be set and ε is the end condition of the iterative solution of the variational model The VMD algorithm decomposes the collected pressure wave signal to obtain the time domain diagram of each BIMF. Then, it gets the spectrum diagram of each BIMF through an FTT, as shown in Figures 9 and 10. The useful low-frequency information is mainly distributed in the low-frequency band, and the noise interference is mainly distributed in the high-frequency band. Therefore, the main frequency component can be selected by referring to the spectrum of each BIMF to reconstruct the signal.

| Signal reconstruction
The number of decomposition layers of the VMD algorithm affects the accuracy of decomposition results, and the noise signal seriously involves extracting helpful information. If the BIMF components obtained from decomposition are directly linearly stacked, the original signal can be accurately F I G U R E 10 Time-frequency diagram of each band-limited intrinsic mode function (BIMF) after pressure wave decomposition when the leakage position is at 2/3 of the annulus.
reconstructed, but this reconstruction cannot achieve the effect of noise removal. Therefore, to solve the problem of the optimal decomposition level of the original signal and the selection of BIMF signal reconstruction, the central frequency method and Hilbert transform are used to analyze the time-frequency and energy characteristics of the signal. Hilbert joint spectrum can combine the three signal characterization scales of time, frequency, and energy, and can more intuitively analyze the signal frequency and power changes with time. VMD algorithm decomposes the original signal into multiple BIMFs and a residual signal, as shown in Equation (24)  For the jth BIMF, its Hilbert transform can be expressed as Equation (25), and the original signal is defined as Equation (26) Integrating H ω t ( , ) and H ω t ( , ) 2 separately over time, the Hilbert marginal spectrum h ω ( ) and Hilbert energy spectrum ES ω ( ) can be obtained as shown in Equations (27) and (28) The BIMF obtained from the VMD decomposition of the pressure wave signal is subjected to the Hilbert transform. Figure 11 shows the frequency and Hilbert joint spectrum of each BIMF when the leakage position is at 1/2 and 2/3 of the wellbore annulus. The abscissa and ordinate on the left side of the Hilbert joint spectrum represent the change of the frequency amplitude of each modal component with time, and the abscissa and ordinate on the right side represent the change of the energy scale of the signal with time. The time-frequency diagram and Hilbert joint spectrum of each BIMF show that the pressure wave signal waveform and the central frequency component are concentrated in the same period. VMD algorithm decomposes the pressure wave signal into low and high-frequency signals. From the joint spectrum of Figure 11, only BIMF1's energy and frequency are concentrated in the same period of excitation waveform. Combined with the frequency spectrum of each BIMF in Figure 10, the signal is reconstructed with BIMF1 as the primary frequency component.
As shown in Figure 12, a comparison is made between the reconstructed signal and the original signal. The original signal waveform has too much burr and severe noise, resulting in distortion. After the VMD algorithm denoising, the signal retains most of the features of the original signal. At the same time, the noise contained in it is effectively removed, the sharp local burrs in the signal are also eliminated, and the signal becomes relatively smooth. The experimental results show that the method can effectively remove the noise and has effectiveness and practicality.

| RESULTS AND DISCUSSION
The experimental wellbore is a straight well, and the leakage valves placed at 18 and 24 m downhole are opened, respectively, and pressure sensors collect the pressure wave signal data. As the experimentally acquired signal data contains too much noise, this noise affects the pressure wave frequency domain feature extraction. It needs to be denoised by the VMD algorithm described above. The denoised signal is then passed through the FFT to obtain the pressure wave signal spectrum. The denoised pressure wave signal spectrum shows that the noise is removed, and the dominant frequency is retained. Figure 13 shows the frequency spectrum of the acquired pressure signal, and the pressure wave response spectrum characteristics are apparent. By comparing the frequency response of a leaking system with that of a system without leaks, the pressure wave, due to the F I G U R E 13 Pressure wave response spectrum for wellbore annulus with and without leaks.
F I G U R E 14 Comparison of pressure wave signal spectrum and pressure head function spectrum when the leakage position is at 1/2 of the annulus. presence of leakage, generates additional resonance points and resonance amplitudes in the pressure main frequency interval of 10-120 Hz. Moreover, the energy loss occurs, and the amplitude decreases. Indicates the presence of a leak near the peak of the main pressure amplitude, which contains characteristic information about the location of the leak point. Figures 14 and 15 compare the pressure head function and de-noised pressure wave spectrum. It can be analyzed that when the leakage position is at 1/2 of the annulus, the characteristics of the wellbore annulus system change. With the horizontal axis frequency increase, the resonance amplitude shows a slow attenuation trend in the low-frequency band. When the leakage position is at 2/3 of the annulus, the characteristics of the wellbore annulus system change. With the horizontal axis frequency increase, the resonance amplitude shows a rapid attenuation trend in the low-frequency band. The results show that the spectral amplitude of the pressure wave signal and the simulated spectral amplitude of the pressure head function is consistent in terms of trend characteristics when leaks occur at different locations in the wellbore annulus so that the location of the leaky layer on the wellbore annulus can be located.
For different leakage amounts, when the leakage valve is kept open, the leakage amount is controlled by F I G U R E 15 Comparison of pressure wave signal spectrum and pressure head function spectrum when the leakage position is at 2/3 of the annulus.
F I G U R E 16 Pressure wave spectra of different leakage amounts when leaking at 1/ 2 of the annulus. using a flow meter. The leakage amount is 6% and 10% of the average flow rate for the experiment. The pressure wave spectra under different leakage conditions are plotted on the same graph, as shown in Figures 16  and 17, from which it can be shown that when the leakage amount increases, the amplitude of the pressure wave spectra decays, which indicates that the larger the leakage amount, the greater the energy loss, resulting in a decrease in the amplitude of the pressure wave spectra. For different leakage amounts at the same position, while keeping the operating characteristics unchanged, the trend of the spectrum curves remains the same, with a decrease in the amplitude of the pressure wave spectrum.

| CONCLUSIONS
The method performs leak detection in the wellbore annulus using an excitation impulse response. The transfer function of the annular wellbore system is only related to its characteristics. Once leaks occur in the wellbore annulus, the transfer function of the wellbore annulus system will change accordingly. Under the condition that the operating features of the control valve remain unchanged, the pressure head function is only related to the system transfer function. Thus, leaks can be located based on the frequency domain curve of the impulse pressure response.
The pressure wave signal collected in the experiment contains too much noise. The VMD algorithm is applied to decompose the collected pressure wave signal. The VMD algorithm combines with the Hilbert joint spectrum to extract the main frequency signal to reconstruct the signal for denoising. When the pressure wave contains a lot of noise, it is difficult to identify the peak of the negative pressure wave when the pressure wave signal is analyzed in the time domain, causing some difficulties in locating the leaky layer. When the pressure wave signal is analyzed in the frequency domain, the frequency domain features are apparent, making it easier to identify the location of the leaky layer.
Compared with using instruments, this method also has additional tools, but it does not need to go down into the well. All devices and operations are done at the surface, and the process is simple, the response is rapid and active excitation, the data acquisition and analysis time is fast, and the pressure wave signal waveform is obvious. According to the experimental verification, it can realize the effective positioning of the leakage location for different depths of leakage. In summary, the principle of the method is feasible and has some practical significance for field application.

AUTHOR CONTRIBUTIONS
Wanneng Lei and Yanxian Wu were involved in data curation. Zhongxi Zhu and Beijing Wang were involved in methodology, numerical simulation, and writingoriginal draft. Yanxian Wu provided supervision. Zhongxi Zhu and Beijing Wang contributed to writing-review and editing. All authors have read and agreed to the published version of the manuscript.

ACKNOWLEDGMENTS
The authors wish to acknowledge the China Scholarship Council and National Engineering Laboratory of Petroleum Drilling Technology, Leak Resistance and Sealing Technology Research Department. This research was F I G U R E 17 Pressure wave spectra of different leakage amounts when leaking at 2/ 3 of the annulus. funded by the Hubei Provincial University Outstanding Young and Middle-aged Science and Technology Innovation Team Plan (Grant No. T201804)