Underwater ultrasonic Doppler real-time analysis method and device based on digital transducer

Abstract

The invention provides an underwater ultrasonic Doppler real-time analysis method and device based on a digital transducer, and relates to the technical field of data analys.The method comprises the steps that 1, multiple beams of directional ultrasonic waves are emitted into drilling fluid through a digital transducer array on the outer wall of a marine riser, echo signals reflected by gas-liquid two-phase flow are received, and the digital transducer array is used for receiving the echo signals; performing real-time digital conversion on the echo signal to generate an original digital signal containing phase and amplitude information; 2, Doppler frequency shift feature extraction is conducted on the original digital signals, and the real-time movement speed and direction of bubbles in the drilling fluid are determined by calculating the frequency deviation of reflected ultrasonic waves and transmitted ultrasonic waves. Through the digital transducer array and the ultrasonic Doppler technology, accurate measurement of motion parameters of bubbles in drilling fluid, real-time calculation of local gas content and early warning and grading response of gas cut are achieved, and the safety and the intelligent level of drilling operation are improved.

Description

Technical Field

[0001] The present invention relates to the technical field of data analysis, and particularly to an underwater ultrasonic Doppler real-time analysis method and device based on digital transducers. Background Art

[0002] In deep-water drilling operations, early monitoring of gas invasion is a key link in preventing blowout accidents. Currently, mainstream monitoring technologies are mostly based on the detection principle of the macroscopic volume change of drilling fluid, mainly including the inlet and outlet flow difference method and the drilling fluid pit liquid level monitoring method. The former judges the gas invasion situation by comparing the difference between the mud pumping volume and the return volume, and the latter calculates the volume of the invaded gas based on the change in the height of the drilling fluid pit liquid level.

[0003] However, in the deep-water riser environment, these methods face a series of technical challenges. In the initial stage of gas invasion, the invaded gas forms millimeter-scale microbubble clusters under high-pressure conditions, which has a relatively small impact on the overall volume of the drilling fluid. Limited by sensitivity, traditional monitoring methods are difficult to capture such microscopic changes and often need the gas invasion to develop to a certain extent before an alarm can be triggered, thus missing some early disposal opportunities.

[0004] In addition, during the ascent of gas in the riser, it will experience complex flow pattern conversions such as bubbly flow, slug flow, and plug flow. The coalescence and breakup of bubbles make it difficult to predict the slip velocity. Traditional methods mostly rely on fixed slip models to calculate the gas migration velocity, and in practical applications, it is easy to cause deviation in the judgment of the gas invasion position due to flow pattern changes. Summary of the Invention

[0005] The technical problem to be solved by the present invention is to provide an underwater ultrasonic Doppler real-time analysis method and device based on digital transducers to achieve early identification, precise positioning, and intelligent prevention and control of drilling fluid gas invasion.

[0006] To solve the above technical problems, the technical solution of the present invention is as follows:

[0007] In the first aspect, an underwater ultrasonic Doppler real-time analysis method based on digital transducers, the method includes:

[0008] Step 1, transmitting multiple beams of directional ultrasonic waves into the drilling fluid through a digital transducer array on the outer wall of the riser, receiving the echo signals reflected by the gas-liquid two-phase flow, and performing real-time digital conversion on the echo signals to generate original digital signals containing phase and amplitude information;

[0009] Step 2, extracting the Doppler frequency shift characteristics of the original digital signals, and determining the real-time movement speed and direction of the bubbles in the drilling fluid by calculating the frequency deviation between the reflected ultrasonic waves and the transmitted ultrasonic waves;

[0010] Step 3: Based on the flow field characteristics in the riser, establish a three-dimensional dynamic spatial grid, spatially discretize and map the real-time movement speed and direction of the bubbles according to the grid positions, and generate a set of bubble movement vectors for each grid cell;

[0011] Step 4: Analyze the gradient change of the set of bubble movement vectors in adjacent grid cells, generate a grid characteristic difference value based on the sudden change of the vector direction and the rate difference, and dynamically correct the set of bubble movement vectors based on the grid characteristic difference value, and output the optimized bubble movement parameters;

[0012] Step 5: Based on the optimized bubble movement parameters, combined with the analysis of the acoustic attenuation characteristics on the ultrasonic propagation path, calculate the local gas holdup distribution in each grid cell in real time;

[0013] Step 6: Conduct a time-series dynamic tracking of the local gas holdup distribution, identify the starting position, diffusion trend and acceleration of the abnormal change of the gas holdup, generate an early warning signal for gas invasion, and perform real-time pattern matching of the warning signal with a preset threshold value to automatically trigger the well control instruction at different levels.

[0014] Furthermore, transmit multiple beams of directional ultrasonic waves into the drilling fluid through a digital transducer array on the outer wall of the riser, receive the echo signals reflected by the gas-liquid two-phase flow, and perform real-time digital conversion on the echo signals to generate an original digital signal containing phase and amplitude information, including:

[0015] Capture the ultrasonic echo analog signals at different axial depth positions through a transducer array arranged annularly on the outer wall of the riser;

[0016] Detect the spectral characteristics of the drilling fluid background noise, remove the components matching the noise spectrum from the ultrasonic echo analog signals, and generate a primary purification signal;

[0017] According to the change of the acoustic characteristics of the drilling fluid in the deep water high-pressure environment, adjust the gain amplitude of the primary purification signal in real time to generate a secondary optimized signal with balanced intensity;

[0018] Discretely sample the secondary optimized signal at a preset time interval, record the amplitude value, phase angle and accurate time stamp of each sampling point to form a digital signal sequence;

[0019] Perform phase alignment and amplitude superposition on the multi-channel digital signals with overlapping time windows and spatial coordinates to generate an original digital signal with three-dimensional position tags.

[0020] Furthermore, extract the Doppler frequency shift characteristics of the original digital signal, and determine the real-time movement speed and direction of the bubbles in the drilling fluid by calculating the frequency deviation between the reflected ultrasonic wave and the transmitted ultrasonic wave, including:

[0021] The original digital signal with three-dimensional position tags is segmented into continuous time windows at a fixed duration. The ultrasonic frequency offset and the phase change trajectory are analyzed within each window, and a time-frequency feature matrix is output.

[0022] Based on the time-frequency feature matrix, signal regions with continuously similar frequency shift amplitudes and consistent phase evolution directions are identified and segmented into independent bubble acoustic fingerprint clusters.

[0023] For each bubble acoustic fingerprint cluster, the difference in frequency shift between adjacent time windows is calculated along the ultrasonic emission direction to generate the axial movement speed component of the bubble. Based on the phase difference between adjacent receiving units in the bubble acoustic fingerprint cluster, the tangential offset angle component of the bubble group is calculated.

[0024] Furthermore, based on the time-frequency feature matrix, signal regions with continuously similar frequency shift amplitudes and consistent phase evolution directions are identified and segmented into independent bubble acoustic fingerprint clusters, including:

[0025] In the time-frequency feature matrix, the frequency shift amplitude sequence of the same spatial coordinate point within three consecutive time windows is extracted, and the absolute difference in frequency shift between adjacent windows is calculated. If the absolute differences are all less than the dynamic fluctuation threshold, the coordinate point is marked as a frequency shift stable point.

[0026] For the cluster of frequency shift stable points, the phase change trajectory is analyzed, and the combination of phase difference signs between adjacent receiving units within the same time window is calculated. When the combination of phase difference signs meets the preset direction consistency condition, it is determined that the phase evolution direction is consistent.

[0027] The cluster of frequency shift stable points that meet the condition of consistent phase evolution direction is clustered into independent bubble acoustic fingerprint clusters according to the principle of spatial continuity.

[0028] Furthermore, a three-dimensional dynamic space grid is established based on the flow field characteristics in the riser. The real-time movement speed and direction of the bubbles are spatially discretely mapped according to the grid positions to generate a set of bubble movement vectors for each grid cell, including:

[0029] The axial movement speed component of the bubble and the tangential offset angle component are fused to generate a three-dimensional bubble movement vector with spatial coordinate tags.

[0030] According to the inner diameter size of the riser and the characteristics of the gas-liquid two-phase flow, the monitoring area is divided into annular grid cells that are stratified axially and sectorized circumferentially, and the three-dimensional bubble movement vectors are assigned to the corresponding annular grid cells according to the spatial coordinate tags.

[0031] The time stamps of the vector data within each grid cell that lasts for a preset duration are verified to generate the current valid set of bubble movement vectors.

[0032] Further, perform gradient change analysis on the bubble motion vector sets of adjacent grid cells, generate grid feature difference values based on the abrupt change of vector direction and rate difference, and dynamically correct the bubble motion vector sets based on the grid feature difference values, and output optimized bubble motion parameters, including:

[0033] For the current valid bubble motion vector set, calculate the maximum deviation value and rate change rate of the motion direction angle between each grid cell and its adjacent cells;

[0034] Generate grid feature difference values based on the maximum deviation value and rate change rate;

[0035] Compare the grid feature difference values with a preset threshold. When the grid feature difference value is greater than the preset threshold, generate a strong correction factor; when the grid feature difference value is less than the preset threshold, generate a weak correction factor;

[0036] Based on the strong correction factor and weak correction factor, perform neighborhood smoothing on the motion vectors of the current grid cell, and output optimized bubble motion parameters.

[0037] Further, based on the optimized bubble motion parameters, in combination with the analysis of the acoustic attenuation characteristics on the ultrasonic propagation path, calculate the local gas holdup distribution in each grid cell in real time, including:

[0038] According to the optimized bubble motion parameters, extract the average motion rate of bubbles and the number of bubbles per unit volume in each grid cell, and calculate the sound intensity attenuation rate before and after the ultrasonic wave penetrates the current grid cell;

[0039] Calculate the additional acoustic attenuation coefficient caused by the bubble group according to the sound intensity attenuation rate and the number of bubbles;

[0040] Determine the gas volume fraction through the correlation relationship between the additional acoustic attenuation coefficient and the equivalent radius of the bubble;

[0041] Real-time collect the temperature, pressure and solid phase concentration of the drilling fluid, compensate the gas volume fraction, and output the local gas holdup distribution value with three-dimensional grid coordinates.

[0042] Further, perform time-series dynamic tracking on the local gas holdup distribution, identify the starting position, diffusion trend and acceleration of abnormal gas holdup changes, generate early gas invasion warning signals, and perform real-time pattern matching of the warning signals with a preset threshold to automatically trigger well control instructions at different levels, including:

[0043] Continuously cache the local gas holdup distribution values of the current and the previous two monitoring periods, and construct a gas holdup spatio-temporal distribution matrix, the matrix dimension includes grid coordinates, gas holdup values, and timestamps;

[0044] In the spatiotemporal evolution matrix, the regions where the gas fraction increment of the same grid cell in adjacent periods is greater than or equal to the first threshold are detected and marked as abnormal source cells;

[0045] Taking the abnormal source unit as the center, the gas content gradient change direction of adjacent grids in the axial or radial direction is analyzed, and the diffusion vector of the abnormal area is predicted by combining the bubble motion vector set;

[0046] Calculate the acceleration of the gas content change in the abnormal area. When the acceleration is greater than or equal to the second threshold, a first-level warning signal is generated. When the diffusion vector of the abnormal area points toward the wellhead and the cumulative increase in gas content is greater than or equal to the third threshold, a second-level warning signal is generated.

[0047] The early warning signals are ranked and the corresponding well control instructions are automatically triggered. For example, the first-level early warning automatically adjusts the drilling pump displacement to compensate for the annular space pressure, and the second-level early warning triggers the pre-closure of the semi-sealed gate of the blowout preventer group and the pressurization of the throttle manifold.

[0048] Furthermore, with the abnormal source unit as the center, the gas fraction gradient change direction of the adjacent grids in the axial or radial direction is analyzed, and the diffusion vector of the abnormal area is predicted in combination with the bubble motion vector set, including:

[0049] Extract the three-dimensional grid coordinates of the abnormal source unit, and obtain the local gas content distribution values of the upstream grid unit and downstream grid unit adjacent to the abnormal source unit in the axial direction, and the circumferential grid unit adjacent to the radial direction;

[0050] Based on the local gas fraction distribution value, the gas fraction difference between the abnormal source unit and each adjacent unit is calculated to generate gas fraction gradient change direction data;

[0051] According to the average bubble motion velocity vector of the grid where the abnormal source unit is located and the adjacent grids in the bubble motion vector concentration, the gas fraction gradient change direction data is vector superimposed and corrected to generate the abnormal area diffusion vector.

[0052] Secondly, an underwater ultrasonic Doppler real-time analysis device based on a digital transducer includes:

[0053] A signal acquisition unit, configured to transmit ultrasonic waves through a digital transducer array on the outer wall of the riser, receive echo signals, and output an original digital signal containing phase and amplitude information;

[0054] Frequency shift analysis unit, used to extract Doppler frequency shift features from the original digital signal and output the real-time movement speed and direction of the bubble;

[0055] A grid mapping unit is used to establish a three-dimensional dynamic space grid and discretize and map the bubble motion parameters, and output a bubble motion vector set for each grid unit;

[0056] A vector correction unit is used to perform gradient analysis and dynamic correction on the bubble motion vector sets of adjacent grid cells, and output optimized bubble motion parameters.

[0057] A gas holdup calculation unit is used to calculate the local gas holdup distribution in each grid cell in real time based on the optimized bubble motion parameters and in combination with the analysis of the acoustic attenuation characteristics on the ultrasonic propagation path.

[0058] An early warning execution unit is used to dynamically track the local gas holdup distribution and trigger well control instructions in a hierarchical manner.

[0059] The above solution of the present invention has at least the following beneficial effects:

[0060] Through the extraction of Doppler frequency shift characteristics, the motion speed and direction of bubbles can be calculated in real time, and the three-dimensional motion trajectory of bubbles in the drilling fluid can be accurately captured, solving the problem that it is difficult to dynamically track tiny bubbles by traditional methods. By establishing a three-dimensional dynamic space grid and mapping the bubble motion parameters according to the grid positions, a refined modeling of the flow field in the riser is realized, which can intuitively reflect the bubble distribution and motion characteristics in different regions and provide basic data for the analysis of gas-liquid two-phase flow. Through gradient analysis and dynamic correction, the parameter mutation and noise interference between adjacent grids are eliminated, making the bubble motion vector sets more conform to the actual flow field law and improving the data stability.

[0061] Combined with the acoustic attenuation characteristics and environmental parameter compensation, considering the effects of temperature, pressure, and solid concentration on the gas volume and ultrasonic propagation, the local gas holdup distribution with three-dimensional coordinates is output in real time, significantly improving the accuracy and reliability of gas holdup monitoring. Through the time series analysis of the gas holdup distribution, the starting position, diffusion direction, and acceleration of gas invasion can be quickly identified, capturing the signs of gas invasion earlier than traditional methods and shortening the early warning time. According to the early warning level, well control instructions are automatically triggered (such as adjusting the pump displacement, pre-closing the blowout preventer, etc.), realizing the full-process automation from early warning to control, reducing the delay of manual intervention, improving the well control efficiency, and effectively reducing the risk of accidents such as blowout caused by gas invasion.

[0062] From the bubble motion parameters to the gas holdup distribution and then to the gas invasion early warning, a complete data chain is formed, providing real-time and comprehensive downhole fluid state information for drilling engineers and assisting in optimizing the drilling fluid circulation parameters and adjusting the operation plan. BRIEF DESCRIPTION OF THE DRAWINGS

[0063] Figure 1 is a flow schematic diagram of an underwater ultrasonic Doppler real-time analysis method based on a digital transducer provided by an embodiment of the present invention.

[0064] Figure 2 is a schematic diagram of an underwater ultrasonic Doppler real-time analysis device based on a digital transducer provided by an embodiment of the present invention. DETAILED DESCRIPTION OF THE INVENTION

[0065] Exemplary embodiments of the present disclosure will be described in more detail below with reference to the accompanying drawings. Although the exemplary embodiments of the present disclosure are shown in the drawings, it should be understood that the present disclosure can be implemented in various forms and should not be limited by the embodiments set forth herein. On the contrary, these embodiments are provided so that the present disclosure can be more thoroughly understood and the scope of the present disclosure can be fully conveyed to those skilled in the art.

[0066] As Figure 1 shown, an embodiment of the present invention provides a real-time underwater ultrasonic Doppler analysis method based on digital transducers, and the method includes the following steps:

[0067] Step 1, transmitting multiple beams of directional ultrasonic waves into the drilling fluid through a digital transducer array on the outer wall of the riser, receiving the echo signals reflected by the gas-liquid two-phase flow, and performing real-time digital conversion on the echo signals to generate original digital signals containing phase and amplitude information;

[0068] Step 2, extracting the Doppler frequency shift characteristics of the original digital signals, and determining the real-time movement speed and direction of the bubbles in the drilling fluid by calculating the frequency deviation between the reflected ultrasonic waves and the transmitted ultrasonic waves;

[0069] Step 3, establishing a three-dimensional dynamic space grid based on the flow field characteristics inside the riser, spatially discretizing and mapping the real-time movement speed and direction of the bubbles according to the grid positions, and generating a set of bubble movement vectors for each grid cell;

[0070] Step 4, performing gradient change analysis on the set of bubble movement vectors of adjacent grid cells, generating grid characteristic difference values according to the sudden change of vector direction and rate difference, and dynamically correcting the set of bubble movement vectors based on the grid characteristic difference values to output optimized bubble movement parameters;

[0071] Step 5, based on the optimized bubble movement parameters, combined with the analysis of the acoustic attenuation characteristics on the ultrasonic wave propagation path, calculating the local gas holdup distribution in each grid cell in real time;

[0072] Step 6, performing time-series dynamic tracking on the local gas holdup distribution, identifying the starting position, diffusion trend and acceleration of abnormal gas holdup changes, generating early gas intrusion warning signals, and performing real-time pattern matching between the warning signals and preset thresholds to automatically trigger well control instructions at different levels.

[0073] In the embodiments of the present invention, a digital transducer array is combined with multi-beam directional ultrasonic emission to achieve three-dimensional spatial signal acquisition of the drilling fluid in the riser, significantly improving the data coverage compared to single-point sensors. The real-time digital conversion and noise rejection mechanism of the echo signal avoid electromagnetic interference and quantization errors in analog signal transmission, effectively improving the signal-to-noise ratio and accuracy of the original data. Through time-frequency feature matrix analysis and bubble acoustic fingerprint cluster segmentation technology, multi-dimensional information such as phase and frequency in the echo signal is deeply mined, and bubble groups in different motion states are accurately separated. Compared with traditional simple frequency shift calculation, it can more accurately determine the real-time motion speed and direction of bubbles, especially suitable for measuring the dynamic parameters of micro-bubbles under complex flow patterns. Based on the flow field characteristics in the riser, a three-dimensional dynamic space grid is constructed, and the bubble motion parameters are spatially discretely mapped to achieve refined modeling of the gas-liquid two-phase flow. This method breaks through the limitations of traditional macroscopic parameter monitoring, can visually present the motion states of bubbles at different positions, and provides a quantitative basis for local flow field analysis.

[0074] Through gradient change analysis and dynamic correction strategies, the bubble motion vectors of adjacent grid cells are optimized, effectively suppressing data mutation noise and flow field interference factors. This mechanism can adapt to the flow pattern conversion problem during gas ascent, improving the stability and reliability of bubble motion parameters and reducing monitoring errors caused by complex flow patterns. Combining bubble motion parameters with ultrasonic acoustic attenuation characteristics and comprehensively considering environmental factors such as the temperature and pressure of the drilling fluid, the real-time calculation of the local gas holdup of each grid cell is achieved. Compared with traditional gas holdup estimation methods that rely on overall volume changes, this technology can capture micro gas invasion events and improve the sensitivity of early gas invasion monitoring. Through time series dynamic tracking and hierarchical threshold matching, the starting position of abnormal gas holdup changes is located, the diffusion trend is predicted, and the acceleration is analyzed. Compared with existing technologies, this early warning mechanism significantly shortens the response time from signal acquisition to early warning output, and can automatically trigger well control commands according to the severity of gas invasion, providing more timely and accurate safety protection for deepwater drilling operations.

[0075] In a preferred embodiment of the present invention, in step 1 above, multi-beam directional ultrasonic waves are emitted into the drilling fluid through a digital transducer array on the outer wall of the riser, and the echo signals reflected by the gas-liquid two-phase flow are received. The echo signals are subjected to real-time digital conversion to generate original digital signals containing phase and amplitude information, which may include:

[0076] Step 100, capturing ultrasonic echo analog signals at different axial depth positions through a transducer array arranged annularly on the outer wall of the riser;

[0077] Step 101, detecting the spectral characteristics of the drilling fluid background noise, removing the components matching the noise spectrum from the ultrasonic echo analog signals, and generating primary purified signals;

[0078] Step 102: According to the change of the acoustic characteristics of the drilling fluid in the deep - water high - pressure environment, adjust the gain amplitude of the primary purification signal in real time to generate a secondary optimized signal with balanced intensity.

[0079] Step 103: Discretely sample the secondary optimized signal at preset time intervals, record the amplitude values, phase angles, and precise timestamps of each sampling point to form a digital signal sequence.

[0080] Step 104: Align the phases and superimpose the amplitudes of the multi - channel digital signals with coincident time windows and spatial coordinates to generate an original digital signal with three - dimensional position tags.

[0081] In the embodiment of the present invention, a circularly arranged digital transducer array is installed on the outer wall of the riser. Assume the array consists of 32 transducers, evenly distributed on the circumference. Each transducer has an independent transmitting and receiving circuit, which can precisely control the time of ultrasonic wave transmission and reception. When transmitting ultrasonic waves, first set the transmission frequency to 500 kHz, and each transducer emits an ultrasonic beam into the drilling fluid at a specific tilt angle. For example, the transducer at the top emits at a downward angle of 30°, and the bottom transducer emits at an upward angle of 30° to ensure that the ultrasonic waves can cover different depth areas in the riser.

[0082] After the transmission is completed, the transducer immediately switches to the receiving mode. According to the propagation speed of ultrasonic waves in the drilling fluid, set the receiving time window. For example, if it is necessary to detect the echo at a distance of 0.5 meters from the transducer, since the ultrasonic wave travels back and forth and the propagation distance is 1 meter, calculate the receiving time. Each transducer captures the ultrasonic echo analog signal within its respective time window, and these signals are recorded by the acquisition device in the form of electrical signals to form the original analog signal data.

[0083] Step 101: The collected analog signals are mixed with various noises, such as mechanical noises generated by the vibration of drilling equipment, thermal noises of circuit components, etc. First, segment the collected signals, and set the duration of each segment to 10 milliseconds. Then, use a spectrum analyzer to perform spectrum analysis on each segment of the signal to obtain its frequency - amplitude distribution. Through the spectrum analysis of a large number of signal segments, statistically analyze the main frequency components of the drilling fluid background noise. For example, it is found that the noise has a higher amplitude in the frequency bands of 200 - 300 Hz and 800 - 1000 Hz. Accordingly, construct a noise spectrum template and record the average amplitudes in these frequency bands.

[0084] When removing noise, the original analog signal is subjected to spectral analysis again. For the parts with high matching degrees of frequency and amplitude in the noise spectral template, suppression processing is performed. The specific judgment method is as follows: If the signal amplitude within a certain frequency band exceeds 1.2 times the amplitude of the corresponding frequency band in the noise template, the signal in this frequency band is retained; otherwise, the amplitude of the signal in this frequency band is reduced to 20% of the original value to eliminate most of the background noise and obtain a primary purified signal.

[0085] Step 102, in a deep - water environment, the properties of drilling fluid such as density and viscosity change with depth, affecting the propagation of ultrasonic waves. Through the pressure sensors and temperature sensors installed on the riser, the pressure and temperature data of the drilling fluid are obtained in real - time. Assume that a relationship table between the acoustic properties of the drilling fluid and signal attenuation is established. For example, when the pressure is 10 MPa and the temperature is 15 °C, the ultrasonic attenuation coefficient is 0.5 dB / m; when the pressure becomes 15 MPa and the temperature is 20 °C, the attenuation coefficient becomes 0.8 dB / m. According to the real - time collected pressure and temperature data, query the relationship table to determine the current ultrasonic attenuation coefficient. Compare it with the preset standard attenuation coefficient (assumed to be 0.3 dB / m) and calculate the gain adjustment amplitude required. If the current attenuation coefficient is 0.8 dB / m, it indicates that the signal attenuation is serious and the gain needs to be increased. Calculate the gain adjustment ratio as 0.8÷0.3≈2.67, that is, amplify the amplitude of the primary purified signal by 2.67 times through the gain amplifier to obtain a secondary optimized signal with balanced intensity.

[0086] Step 103, connect the secondary optimized signal to a high - speed sampling device, and set the sampling frequency to 10 MHz, that is, sample once every 0.1 microseconds. At each sampling moment, the acquisition device measures the instantaneous voltage value of the signal, and this voltage value corresponds to the amplitude of the signal. To obtain the signal phase angle, compare the sampled signal with a reference signal that is of the same frequency and in - phase as the transmitted ultrasonic wave. When the time difference between the rising edge of the sampled signal and the rising edge of the reference signal is 0.05 microseconds (corresponding to a 180° phase difference), determine that the phase angle of the sampled signal is 180°; if the time difference is 0.025 microseconds, the phase angle is 90°, and so on, calculate the phase angle of each sampling point. At the same time, record the exact time of each sampling point, accurate to the nanosecond level. Organize the amplitude, phase angle, and timestamp information of each sampling point into a data record, and arrange them in sequence according to the sampling order to form a digital signal sequence. For example, the first sampling point is recorded as (amplitude: 2.3 V, phase angle: 60°, timestamp: 00:00:00.000000123), the second sampling point is recorded as (amplitude: 2.1 V, phase angle: 75°, timestamp: 00:00:00.000000223), and so on.

[0087] Step 104: Since 32 transducers collect 32 signals, for multi-channel signals within the same time window (assumed to be 1 millisecond) and with coincident spatial coordinates (i.e., the same depth position), the first signal is used as a reference. Calculate the phase differences between the other 31 signals and the reference signal within this time window. For example, if the rising edge of the second signal is 0.01 microseconds later than the reference signal, it indicates a phase difference of 36°. Through digital signal processing algorithms, perform time delay processing on the second signal, delaying it by 0.01 microseconds to align its phase with the reference signal. Repeat this operation for all 31 signals to complete phase alignment. After phase alignment, perform amplitude superposition. Add the amplitude values of the 32 signals at the same sampling point to obtain the superimposed amplitude. For example, at a certain sampling point, the amplitudes of the 32 signals are 2.3V, 2.1V, 2.0V..., adding these values gives a superimposed amplitude of 65.4V. Combine the position of the transducer in the circular array and the signal acquisition depth to assign a three-dimensional position tag to the superimposed signal. Assume that the 16th transducer is located on the right side of the circumference, with its planar coordinates being (1, 0) (taking the center of the riser as the origin) and the acquisition depth being 10 meters, then the three-dimensional position tag of this signal is (1, 0, 10). Finally, generate the original digital signal with three-dimensional position tags.

[0088] The transducer array arranged in a ring can capture echo signals at different axial depths, enabling three-dimensional spatial monitoring of the drilling fluid in the riser. Compared with the single-point monitoring method, the data obtained is more comprehensive, avoiding monitoring blind spots. Through the noise rejection and gain adjustment steps, background noise interference is effectively removed, adapting to the changes in the acoustic characteristics of the drilling fluid in the deep-water environment, ensuring the stability of the echo signal intensity, and improving the signal-to-noise ratio and reliability of the signal. The discrete sampling process accurately records the amplitude, phase, and time information of the signal. The phase alignment and amplitude superposition operations of the multi-channel signals enhance the signal intensity, and at the same time assign a three-dimensional position tag to the signal, making the original digital signal not only contain rich physical information but also have accurate spatial positioning, improving the accuracy and effectiveness of signal processing.

[0089] In a preferred embodiment of the present invention, in step 2 above, for the extraction of the Doppler frequency shift characteristics of the original digital signal, by calculating the frequency deviation between the reflected ultrasonic wave and the transmitted ultrasonic wave to determine the real-time movement speed and direction of the bubbles in the drilling fluid, it may include:

[0090] Step 200: Divide the original digital signal with three-dimensional position tags into continuous time windows at a fixed duration, and analyze the ultrasonic frequency offset and phase change trajectory within each window, and output a time-frequency feature matrix;

[0091] Step 201: Based on the time-frequency feature matrix, identify signal regions with continuously similar frequency shift amplitudes and the same phase evolution direction, and divide them into independent bubble acoustic fingerprint clusters, specifically including:

[0092] Step 2010, in the time-frequency feature matrix, extract the frequency shift amplitude sequences of the same spatial coordinate point within three consecutive time windows, calculate the absolute difference of the frequency shift amounts between adjacent windows. If the absolute differences are all less than the dynamic fluctuation threshold, mark this coordinate point as a frequency shift stable point;

[0093] Step 2011, cluster the frequency shift stable points, analyze the phase change trajectory, calculate the phase difference symbol combinations between adjacent receiving units within the same time window. When the phase difference symbol combinations meet the preset direction consistency condition, it is determined that the phase evolution directions are consistent;

[0094] Step 2012, cluster the frequency shift stable points that meet the consistent phase evolution directions into independent bubble acoustic fingerprint clusters according to the principle of spatial continuity;

[0095] Step 202, for each bubble acoustic fingerprint cluster, calculate the difference in the frequency shift amount between adjacent time windows along the ultrasonic wave emission direction to generate the axial motion rate component of the bubble; calculate the tangential offset angle component of the bubble group according to the phase difference of the bubble acoustic fingerprint cluster between adjacent receiving units.

[0096] In the embodiment of the present invention, the original digital signal with three-dimensional position tags (circumferential angle, radial distance, axial depth) is segmented into consecutive time windows with a fixed duration (such as 200 milliseconds), and adjacent windows can overlap by 50 milliseconds to ensure signal continuity. For example, the first window covers 0 - 200 milliseconds, the second window covers 100 - 300 milliseconds, and so on. Within each time window, compare the signal frequency with the transmitted ultrasonic wave frequency (such as 800 kHz). Determine the frequency shift by detecting the periodic change of the signal waveform. If the period of the transmitted signal is 1.25 microseconds (corresponding to 800 kHz), and the period of a certain segment of the signal within the window is 1.3 microseconds, then the frequency drops to about 769 kHz, and the shift amount is -31 Hz. Record the frequency shift values at each moment within the window point by point. Based on the signal phase at the start moment of the window (such as 0°), measure the time delay of the subsequent signal waveforms relative to the reference, and convert it into a phase angle. For example, if the signal peak at a certain moment is delayed by 0.05 microseconds, the corresponding phase change is (0.05÷1.25)×360° = 14.4°, and record it as the phase change trajectory. Arrange the data such as the frequency shift amount, the starting phase value, and the phase change amount of each window in the form of "row - time window, column - characteristic parameter" to generate the time-frequency feature matrix.

[0097] Step 2010, in the time-frequency feature matrix, for a certain spatial coordinate point (such as 60° circumferentially, 0.5 meters radially, 10 meters axially), extract the frequency shift amplitude sequences of three consecutive time windows, assumed to be 18 Hz, 17 Hz, and 19 Hz. Calculate the absolute differences of the frequency shift amounts between adjacent windows, that is, |17 - 18| = 1 Hz, |19 - 17| = 2 Hz.

[0098] The preset dynamic fluctuation threshold is 3 Hz (based on historical data statistics). If the differences are all less than the threshold (e.g., 1 Hz and 2 Hz are less than 3 Hz), then mark this point as a frequency shift stable point.

[0099] Step 2011, within the same time window, take the phase values of adjacent receiving units (such as transducers A and B). Assume that A is 30° and B is 45°, the phase difference is 15°, and the sign is positive (B is greater than A). Continuously observe 5 time windows. If the signs of the phase differences between adjacent units are always positive (such as +, +, +, +, +) or show a periodic pattern (such as +, +, -, -, +), then it is determined that the phase evolution directions are consistent; if the signs are irregular (such as +, -, +, -), then it is determined that they are inconsistent.

[0100] Step 2012, in the three-dimensional space of the time-frequency feature matrix, check whether the frequency shift stable points are adjacent in the circumferential direction (such as the angular difference ≤ 30°) and have similar axial depths (such as the difference ≤ 0.5 m). Cluster the points that satisfy spatial adjacency and consistent phase evolution directions into independent bubble acoustic fingerprint clusters. For example, the stable point cluster within 50° - 70° in the circumference and 9 - 10 m in the axial direction, if the phase signs are consistent, then they are clustered into one cluster.

[0101] Step 202, along the ultrasonic wave emission direction, take the frequency shift amounts of a certain fingerprint cluster in two adjacent windows (such as 22 Hz in the n window and 27 Hz in the n + 1 window), and the difference is 5 Hz. According to the preset corresponding relationship (0 - 10 Hz corresponds to 0.01 m / s / Hz), calculate the axial velocity as 5 × 0.01 = 0.05 m / s. Take the phase difference of the fingerprint cluster between adjacent receiving units (such as 75°), and combine it with the unit interval angle (such as 20°), and convert it into the tangential offset angle proportionally. Considering the propagation path correction (such as fluid refraction), finally determine the offset angle as 30°.

[0102] By integrating frequency and phase information through the time-frequency feature matrix, breaking through the limitations of single frequency shift analysis, it can capture the subtle changes in bubble motion and improve the signal analysis accuracy under complex flow patterns. Based on the dual screening of frequency shift stability and phase consistency, it can effectively distinguish bubble groups in different motion states, avoid signal aliasing interference, and is especially suitable for micro-bubbles and multiphase flow environments. The integrated calculation of axial velocity and tangential angle realizes the spatial vector description of bubble motion, provides key data for the three-dimensional modeling of gas-liquid two-phase flow, and supports the analysis of local flow field characteristics (such as bubble rising trajectory, rotation trend). Through the dynamic threshold and spatial clustering algorithm, it reduces noise interference, ensures the real-time performance and reliability of motion parameters, and provides accurate basis for gas invasion warning, drilling fluid circulation optimization, etc.

[0103] In a preferred embodiment of the present invention, in step 3, a three-dimensional dynamic spatial grid is established based on the flow field characteristics in the riser, and the real-time movement speed and direction of the bubbles are discretely mapped spatially according to the grid positions to generate a set of bubble movement vectors for each grid cell, which may include:

[0104] Step 300, fuse the axial movement rate component and the tangential offset angle component of the bubbles to generate a three-dimensional bubble movement vector with spatial coordinate tags;

[0105] Step 301, according to the inner diameter size of the riser and the characteristics of the gas-liquid two-phase flow, divide the monitoring area into annular grid cells that are stratified axially and sectorized circumferentially, and distribute the three-dimensional bubble movement vectors to the corresponding annular grid cells according to the spatial coordinate tags;

[0106] Step 302, perform timestamp verification on the vector data that persists for a preset duration in each grid cell to generate the current valid set of bubble movement vectors.

[0107] In the embodiment of the present invention, the axial movement rate component and the tangential offset angle component corresponding to each independent bubble acoustic fingerprint cluster are obtained. For example, the axial movement rate component of a certain bubble acoustic fingerprint cluster is 0.08 m / s, and the tangential offset angle component is 40°. At the same time, the spatial coordinates where this fingerprint cluster is located are 90° in the circumferential direction of the riser, a radial distance of 0.6 meters, and an axial depth of 12 meters. The axial movement rate component, the tangential offset angle component, and the spatial coordinates are fused. The axial movement rate component reflects the movement speed of the bubble along the ultrasonic emission direction, the tangential offset angle component represents the offset trend of the bubble in the circumferential direction, and the spatial coordinates clarify the position where the bubble is located. Integrating these three parts of information together forms a three-dimensional bubble movement vector that completely describes the movement state of the bubble. This vector contains the speed magnitude, direction, and specific spatial position information of the bubble movement.

[0108] Step 301 divides the monitoring area based on the actual inner diameter of the riser (assuming a riser inner diameter of 1.5 meters) and the characteristics of the gas-liquid two-phase flow. In the axial direction, considering monitoring accuracy and computational efficiency, stratification is performed at regular intervals (e.g., every 0.3 meters). This divides the axial depth into multiple levels, for example, 0-0.3 meters as the first level, 0.3-0.6 meters as the second level, and so on. In the circumferential direction, the circumference is divided into several sectors based on the required accuracy. For example, to more precisely monitor bubble motion in the circumferential direction, the circumference is divided into 24 sectors, each with an angle of 360° ÷ 24 = 15°. This axial and circumferential stratification creates an annular grid of axially layered and circumferentially sectored cells. For each generated 3D bubble motion vector, its associated grid cell is determined based on its spatial coordinate label. For example, if the spatial coordinates of a 3D bubble motion vector are 120° in the circumferential direction, 0.8 meters in the radial direction, and 3.5 meters in the axial depth, it will be assigned to the annular grid cell corresponding to the 3.3-3.6 meter axial layer and the 105-120° circumferential sector. This method achieves a discretized mapping of the bubble motion vector to the 3D spatial grid.

[0109] Step 302: To ensure that the bubble motion vector set can reflect the current bubble motion state in real time, the vector data in each grid unit is timestamp checked. A time duration is preset (assuming it is 2 seconds), and the timestamp of each vector data is checked to determine whether the vector data is more than the preset time duration from the current moment. For example, there are multiple vector data in a grid unit, and the timestamp of one vector shows that its acquisition time is 3 seconds ago, which exceeds the preset 2-second duration. This means that the data can no longer accurately reflect the current bubble motion state, and it is removed from the vector data of the grid unit; the vector data with a timestamp within 2 seconds is retained. After checking and screening the vector data in each grid unit one by one, the current valid bubble motion vector set is finally generated.

[0110] By constructing a three-dimensional dynamic space grid, the monitoring area in the riser is carefully divided. Combining with the spatial discretization mapping of the bubble motion vector, it can visually and accurately present the motion states of bubbles at different positions. Compared with traditional macroscopic monitoring methods, this approach can deeply analyze the characteristics of the local flow field, provide strong support for establishing a high-precision gas-liquid two-phase flow model, and contribute to understanding the complex laws of bubble motion in drilling fluid. The bubble motion vectors are assigned to grid cells, and an effective vector set is generated through timestamp verification, realizing the structured management of data. When conducting data analysis, the bubble motion data within specific grid cells can be quickly located and retrieved, facilitating the comparison and trend analysis of bubble motion in different regions, greatly improving the efficiency of data processing and analysis, and providing timely data support for decision-making in drilling operations. Timestamp verification is performed on the vector data to promptly eliminate stale data, ensuring that the bubble motion vector set reflects the real situation at the current moment. This mechanism effectively avoids analysis errors caused by using outdated data, improves the accuracy and reliability of bubble motion parameter monitoring, and plays an important role in key aspects such as early gas invasion warning and optimization of the drilling fluid circulation system, ensuring the safe and stable progress of drilling operations.

[0111] In a preferred embodiment of the present invention, for step 4, gradient change analysis is performed on the bubble motion vector sets of adjacent grid cells. According to the abrupt change of the vector direction and the rate difference, a grid characteristic difference value is generated, and based on the grid characteristic difference value, the bubble motion vector set is dynamically corrected, and the optimized bubble motion parameters are output, which may include:

[0112] Step 400, for the current effective bubble motion vector set, calculate the maximum deviation value and the rate of change of the motion direction angle between each grid cell and its adjacent cells;

[0113] Step 401, generate a grid characteristic difference value based on the maximum deviation value and the rate of change of the motion direction angle;

[0114] Step 402, compare the grid characteristic difference value with a preset threshold. When the grid characteristic difference value is greater than the preset threshold, a strong correction factor is generated; when the grid characteristic difference value is less than the preset threshold, a weak correction factor is generated;

[0115] Step 403, based on the strong correction factor and the weak correction factor, perform neighborhood smoothing on the motion vectors of the current grid cell, and output the optimized bubble motion parameters.

[0116] In the embodiments of the present invention, for the current valid set of bubble motion vectors, the adjacent relationships of the grid cells in the riser are first clarified. In a three-dimensional dynamic space grid, each grid cell has adjacent cells in the axial, circumferential, and radial directions. For example, for a grid cell in the 3rd axial layer and the 5th circumferential sector, its adjacent cells include the adjacent sectors in the same layer and the corresponding sector cells in the upper and lower layers. For each grid cell, calculate the maximum deviation value of the motion direction angle between it and its adjacent cells.

[0117] Taking a certain grid cell A as an example: The direction angle of the bubble motion vector of cell A is 60°, the direction angle of the adjacent cell B is 80°, the direction angle of cell C is 75°, and the direction angle of cell D is 100°. By comparing the absolute values of the differences in the direction angles between A and B, C, and D, that is, |80° - 60°| = 20°, |75° - 60°| = 15°, |100° - 60°| = 40°, it is obtained that the maximum deviation value of the motion direction angle between cell A and its adjacent cells is 40°.

[0118] When calculating the rate change rate, taking grid cell A as an example again, its bubble motion rate is 0.1 m / s, and the rate of the adjacent cell B is 0.12 m / s. According to the formula (rate of adjacent cell - rate of this cell) ÷ rate of this cell, the rate change rate of cell A relative to cell B can be obtained as (0.12 - 0.1) ÷ 0.1 = 0.2 (i.e., 20%). For each grid cell, calculate the maximum deviation value of the direction angle and the rate change rate with all adjacent cells in the above manner.

[0119] Step 401: Comprehensively process the maximum deviation value of the motion direction angle and the rate change rate calculated for each grid cell to generate a grid characteristic difference value. Here, a weighted fusion method is adopted. According to actual experience and the characteristics of the drilling fluid flow field, different weights are assigned to the direction angle deviation value and the rate change rate. Assume that the weight of the direction angle deviation value is 0.6 and the weight of the rate change rate is 0.4.

[0120] Taking grid cell A as an example: Its maximum deviation value of the direction angle is 40°, and the rate change rate is 20%. First, normalize the direction angle deviation value. Assume that the maximum possible deviation value is 90°, then the normalized direction angle deviation value is 40° ÷ 90° ≈ 0.44. Multiply the normalized direction angle deviation value by the weight 0.6, getting 0.44 × 0.6 = 0.264; multiply the rate change rate by the weight 0.4, getting 0.2 × 0.4 = 0.08. Finally, add the two results, 0.264 + 0.08 = 0.344, and this value is the grid characteristic difference value of grid cell A. Repeat this process for all grid cells to obtain their respective grid characteristic difference values.

[0121] Step 402: Preset a threshold value, which is determined by analyzing a large amount of drilling data and assumed to be 0.2. Compare the grid property difference value of each grid cell with the preset threshold. When the grid property difference value of a grid cell is greater than 0.2, it indicates that the bubble motion characteristics of this grid cell are quite different from those of adjacent cells, and there may be data anomalies or flow field mutations. At this time, a strong correction factor is generated. The value of the strong correction factor is determined according to the degree of difference. For example, if the difference value is 0.344 (such as grid cell A), the strong correction factor can be set to 0.8, indicating that a relatively large correction needs to be made to the motion vector of this grid cell.

[0122] When the grid property difference value is less than 0.2, it shows that the bubble motion characteristics of this grid cell are relatively close to those of adjacent cells, and the data is relatively stable. A weak correction factor is generated. For example, if the difference value of a certain grid cell is 0.15, the weak correction factor can be set to 0.2, that is, only a relatively small correction needs to be made.

[0123] Step 403: After obtaining the strong and weak correction factors of the grid cell, the motion vector of the current grid cell needs to be smoothed in the neighborhood to optimize the bubble motion parameters. Taking grid cell A as an example, its motion vector includes two key pieces of information: the direction angle and the speed. The specific processing process is as follows:

[0124] Direction angle correction:

[0125] Use the method of neighborhood averaging to adjust the direction angle of grid cell A. First, include the direction angles of cell A and adjacent cells B, C, and D in the calculation range. According to the strong correction factor of cell A (a relatively large value indicates a large correction required), assign a higher weight to the direction angle of cell A itself, and a lower weight to the direction angles of adjacent cells. When operating specifically, superimpose the direction angles of each cell according to the corresponding weights. The direction angle with a higher weight has a greater impact on the final result. Then divide the superimposed value by the sum of all weights. Through this weighted averaging method, the direction angle of cell A approaches the average direction angle of adjacent cells, reducing the direction mutation caused by abnormal data and obtaining the corrected direction angle.

[0126] Speed correction:

[0127] For the speed correction of grid cell A, the weighted averaging strategy is also used. Combine the correction factor of cell A to assign different weights to the speed of cell A itself and the speeds of adjacent cells B, C, and D. When the correction factor is relatively large, the weight of the speed of cell A itself is relatively reduced, and the weights of the speeds of adjacent cells are correspondingly increased. After superimposing the speeds of each cell according to the weights, divide by the sum of all weights, so that the speed value of cell A is closer to the average speed level of adjacent cells, eliminating the speed anomaly caused by local data fluctuations, and thus obtaining the corrected speed.

[0128] By double-correcting the direction angle and rate, the motion vector of grid cell A is effectively smoothed, data mutations and abnormal fluctuations are reduced, and finally optimized bubble motion parameters are output.

[0129] Through gradient change analysis and dynamic correction, the mutation noise and outliers in the bubble motion vector data are effectively suppressed. In the complex flow field of drilling fluid, data fluctuations caused by factors such as bubble aggregation and flow pattern mutation are avoided, and the monitored bubble motion parameters can more truly reflect the actual situation. It can dynamically adjust the correction strategy according to the characteristic differences between grid cells. The optimized bubble motion parameters make the gas-liquid two-phase flow model more accurate, which helps to deeply analyze the local flow field characteristics in the riser. For example, accurately judging the motion trend and aggregation area of bubbles, etc., provides more powerful support for optimizing the design of the drilling fluid circulation system and preventing gas invasion accidents, and improves the safety and efficiency of drilling operations.

[0130] In a preferred embodiment of the present invention, in step 5 above, based on the optimized bubble motion parameters, combined with the analysis of the acoustic attenuation characteristics on the ultrasonic propagation path, the local gas holdup distribution in each grid cell is calculated in real time, which may include:

[0131] Step 500, according to the optimized bubble motion parameters, extract the average motion rate of bubbles in each grid cell and the number of bubbles per unit volume, and calculate the sound intensity attenuation rate before and after the ultrasonic wave penetrates the current grid cell;

[0132] Step 501, calculate the additional sound attenuation coefficient caused by the bubble group according to the sound intensity attenuation rate and the number of bubbles;

[0133] Step 502, determine the gas volume fraction through the correlation between the additional sound attenuation coefficient and the equivalent radius of the bubble;

[0134] Step 503, collect the temperature, pressure and solid concentration of the drilling fluid in real time, compensate the gas volume fraction, and output the local gas holdup distribution value with three-dimensional grid coordinates.

[0135] In the embodiment of the present invention, each grid cell is processed based on the optimized bubble motion parameters. The average motion rate of bubbles is extracted from the parameters, and this value reflects the overall movement speed of the bubbles in the grid; at the same time, the number of bubbles per unit volume is determined, which represents the density of the bubbles in the grid cell. When calculating the sound intensity attenuation rate before and after the ultrasonic wave penetrates the current grid cell, the initial sound intensity value before the ultrasonic wave enters the grid cell needs to be obtained first. This value can be obtained through the set parameters when the transducer emits ultrasonic waves and the previous calibration data. When the ultrasonic wave penetrates the grid cell, the transducer receives the echo signal, and the signal processing device analyzes the intensity of the received echo signal to obtain the sound intensity value after penetration.

[0136] Subtract the sound intensity value after penetration from the initial sound intensity value to obtain the reduction in sound intensity. Then divide the reduction in sound intensity by the initial sound intensity value to calculate the sound intensity attenuation rate before and after the ultrasonic wave penetrates the grid cell. For example, if the initial sound intensity is 100 units and the sound intensity after penetration becomes 80 units, the reduction in sound intensity is 20 units, and the sound intensity attenuation rate is 20 units divided by 100 units, which is 20%. In this way, the sound intensity attenuation rate of each grid cell is calculated.

[0137] Step 501, after obtaining the sound intensity attenuation rate and the number of bubbles per unit volume of each grid cell, enter the calculation process of the additional sound attenuation coefficient. This process needs to be based on a pre-established mathematical relationship model, combined with the actual values of two key parameters, and the result is obtained through multiple steps of analysis and operations.

[0138] First, confirm the two core parameters of each grid cell: the sound intensity attenuation rate and the number of bubbles per unit volume. Take a certain grid cell as an example. Suppose its sound intensity attenuation rate is measured and calculated to be 25%, that is, after the ultrasonic wave penetrates the grid cell, the sound intensity is reduced to 75% of the initial value; the number of bubbles per unit volume is obtained through statistical analysis of the acoustic fingerprint clusters of the bubbles, which is 10,000 bubbles per cubic meter. Since the parameter value ranges of different grid cells may vary greatly, these two parameters need to be standardized. For the sound intensity attenuation rate, divide its value by 100% to convert it into a dimensionless value between 0 and 1. The standardized sound intensity attenuation rate of the above grid cell is 0.25; for the number of bubbles per unit volume (suppose it is 100,000 per cubic meter) for normalization, the standardized number of bubbles per unit volume of this grid cell is 10,000÷100,000 = 0.1.

[0139] The pre-established mathematical relationship model clarifies the relationship between the sound intensity attenuation rate, the number of bubbles per unit volume, and the additional sound attenuation coefficient. In the form of a multivariate function, the standardized sound intensity attenuation rate and the number of bubbles per unit volume are used as independent variables. Substitute the standardized sound intensity attenuation rate of 0.25 and the number of bubbles per unit volume of 0.1 of the above grid cell into the model. Inside the model, according to the set operation rules, these two independent variables are processed. For example, the model first performs a weighted sum of the two independent variables, where the weight of the sound intensity attenuation rate is 0.6 and the weight of the number of bubbles per unit volume is 0.4, that is, perform the operation of 0.25×0.6 + 0.1×0.4 to obtain a preliminary result of 0.15 + 0.04 = 0.19.

[0140] The preliminary calculation result is only an intermediate value and needs to be adjusted according to the correction rules preset in the model. The correction rules consider the combined effects of other influencing factors such as drilling fluid viscosity and ultrasonic frequency on acoustic attenuation. Suppose there is a correction factor related to the drilling fluid viscosity in the model. When the current drilling fluid viscosity is in a certain range, the correction factor is 1.2. Multiply the preliminary calculation result 0.19 by the correction factor 1.2. This final value is the additional acoustic attenuation coefficient caused by the bubble group in this grid cell, which comprehensively reflects the influence degree of the acoustic intensity attenuation rate, the number of bubbles per unit volume, and other environmental factors on the attenuation of ultrasonic waves in this grid cell.

[0141] For all grid cells, the additional acoustic attenuation coefficients are calculated one by one according to the three steps of parameter standardization, substituting into the model for calculation, and result correction. Since the acoustic intensity attenuation rate and the number of bubbles per unit volume of each grid cell are different and may be affected by different environmental factors, the finally obtained additional acoustic attenuation coefficients are also different.

[0142] Step 502, set multiple working condition combinations according to the complex environments that may be encountered in actual drilling. If the drilling fluid density has 3 gradients, which are 1.2 g / cm³, 1.4 g / cm³, and 1.6 g / cm³ respectively; the temperature covers the range of 40°C - 80°C, and a test point is set every 10°C; the pressure ranges from 10 MPa to 30 MPa and is adjusted at intervals of 5 MPa; the bubble concentration forms three states of low concentration (1000 bubbles per cubic meter), medium concentration (5000 bubbles per cubic meter), and high concentration (10000 bubbles per cubic meter) by controlling the gas injection volume. Through permutation and combination, a total of 27 different working condition conditions are constructed.

[0143] Under each working condition, ultrasonic waves are emitted through a digital transducer array installed on the outer wall of the riser. The transducer first emits a directional ultrasonic beam at a frequency of 800 kHz. After the ultrasonic waves penetrate the drilling fluid containing bubbles, the reflected echo signals are received. By comparing the intensity of the transmitted signal with that of the received signal, the sound intensity attenuation rate of the ultrasonic waves penetrating the current drilling fluid area is calculated. Then, combined with the reference value of the sound intensity attenuation without bubbles (previously measured in a pure drilling fluid environment), the additional sound attenuation coefficient caused by the bubble group under this working condition is obtained. For example, under a certain working condition, the measured sound intensity attenuation rate is 30%, and the attenuation rate without bubbles is 5%, then the additional sound attenuation coefficient is 30% - 5% = 25%. The equivalent radius of the bubbles is measured by combining microscopic imaging and high-speed photography. A transparent observation window is set on the flow path of the drilling fluid, and a high-speed camera (frame rate of 1000 frames per second) is used to record the movement process of the bubbles. At the same time, the bubbles are magnified and imaged locally through a microscopic lens (magnification factor of 100 times). After the shooting is completed, the video data is processed frame by frame using image analysis software. The software automatically identifies the bubble contours, approximates the irregular bubble shapes as spheres, calculates the diameter of the equivalent circle by inversely deducing from the contour perimeter, and then takes half of the diameter as the equivalent radius of the bubbles. For example, the measured contour perimeter of a certain bubble is 1.57 mm, the equivalent diameter is calculated to be 0.5 mm through the circle perimeter formula, and the equivalent radius is 0.25 mm. A large number of bubbles (100 bubble samples are selected for each working condition) under each working condition are measured, and the average value is taken as the equivalent radius of the bubbles under this working condition.

[0144] The data of the additional sound attenuation coefficient and the equivalent radius of the bubbles collected under 27 working conditions are sorted into a table form. The parameter combinations of different working conditions (such as drilling fluid density, temperature, pressure, bubble concentration) are listed horizontally in the table, and the corresponding values of the additional sound attenuation coefficient and the equivalent radius of the bubbles are recorded vertically to form a complete data set. The data is initially visualized. Taking the additional sound attenuation coefficient as the abscissa and the equivalent radius of the bubbles as the ordinate, a scatter plot is drawn. By observing the distribution trend of the scatter points, it is found that the data points generally show a linear distribution characteristic, that is, as the additional sound attenuation coefficient increases, the equivalent radius of the bubbles also shows an upward trend, and it is initially judged that there may be a linear relationship between the two. The least squares method is used for linear regression analysis to minimize the sum of the squares of the vertical distances from all data points to this straight line. Specifically, when operating, the data of the additional sound attenuation coefficient and the equivalent radius of the bubbles are imported into professional data analysis software (such as Origin), and a linear regression model is selected for calculation. The software automatically adjusts the slope and intercept of the straight line. After multiple iterative calculations, the equation of the fitted straight line is finally obtained: equivalent radius of the bubbles = 0.5 × additional sound attenuation coefficient + 0.1.

[0145] Taking a certain grid cell as an example, assume that the additional sound attenuation coefficient of this grid cell is calculated to be 0.3 through step 501. Substitute this value into the above correlation relationship, and the equivalent radius = 0.5×0.3 + 0.1 = 0.25 (the unit is assumed to be millimeters). In this way, the equivalent radius of the bubbles in this grid cell is obtained. After knowing the equivalent radius of the bubbles, combined with the number of bubbles per unit volume (assuming the number of bubbles per unit volume in this grid cell is 8000 per cubic meter), calculate the volume of a single bubble according to the sphere volume formula (sphere volume = 4÷3×π×radius³), and then multiply the volume of a single bubble by the number of bubbles per unit volume to obtain the total volume of bubbles per unit volume. Finally, divide the total volume of bubbles per unit volume by the total volume of the grid cell (assuming the volume of the grid cell is 1 cubic meter) to obtain the gas volume fraction, and this value represents the proportion of the gas volume in the total volume within this grid cell.

[0146] Step 503, through sensors installed at different positions of the riser, collect environmental parameters such as the temperature, pressure, and solid concentration of the drilling fluid in real time. The temperature sensor uses a high-precision thermocouple sensor to collect data once per second; the pressure sensor uses a piezoresistive pressure sensor and also maintains a sampling frequency of once per second; the solid concentration is measured by a laser particle size analyzer and updated once every 10 seconds. For example, the data collected at a certain moment is: temperature 60°C, pressure 15 MPa, and solid concentration 12%.

[0147] The pre-established compensation model is constructed based on theories such as thermodynamics and fluid mechanics, and combined with a large amount of actual drilling data. When calculating the local gas holdup distribution, the compensation model needs to comprehensively correct the gas volume fraction considering factors such as temperature, pressure, and solid concentration. The specific analysis process is as follows:

[0148] Substitute the collected temperature T into the compensation model. Based on the ideal gas state equation PV = mRT (where P is the pressure, V is the volume, m is the amount of substance, and R is the gas constant), when the pressure P and the amount of substance m remain unchanged, the change in temperature T will cause a change in the gas volume V.

[0149] For example, it is known that the initial temperature T1 of a certain grid cell is 293 K (corresponding to 20°C), and the gas volume V1 is 1 cubic meter. When the temperature rises to T2 of 333 K (corresponding to 60°C), according to the formula V2 = V1×(T2÷T1), V2 is assumed to be 1.136 cubic meters, and V2 - V1 = 0.136, that is, the gas volume expands by 0.136 cubic meters. This will directly lead to an increase in the gas volume fraction. Assuming the total volume of the original grid cell is 1 cubic meter, only the increase in temperature causes the gas volume fraction to increase by 0.136÷1×100%. Record this increment for subsequent comprehensive calculations.

[0150] The pressure P has a significant impact on the solubility of gas in drilling fluid. The model is calculated based on Henry's law (the solubility of gas is proportional to the pressure). Assuming that the solubility coefficient of a certain gas in drilling fluid is 0.05 (unit: the volume ratio of gas dissolved in unit volume of drilling fluid per megapascal), when the pressure increases from P1 = 10 MPa to P2 = 15 MPa, the volume of gas dissolved more in unit volume of drilling fluid is 0.05×(15-10)=0.25 cubic meters. After this part of the gas is dissolved, the gas phase volume decreases, thus reducing the gas phase volume fraction. If the volume of the grid unit is 1 cubic meter, the gas phase volume fraction is reduced by 0.25÷1×100% due to the pressure change, and this change amount is also recorded for comprehensive calculation.

[0151] The solid phase concentration will affect the propagation characteristics of ultrasonic waves, and thus indirectly affect the gas phase volume fraction. The compensation model incorporates the solid phase concentration - ultrasonic attenuation correction relationship established through a large number of experiments. When the collected solid phase concentration is 12%, according to the correction relationship table or pre-designed calculation rules, it is determined that the additional acoustic attenuation coefficient calculated in step 502 needs to be multiplied by a correction factor of 1.1. After the additional acoustic attenuation coefficient is changed, it is substituted into the correlation relationship between the additional acoustic attenuation coefficient and the equivalent radius of the bubble to recalculate the equivalent radius of the bubble, and then a new gas phase volume fraction is obtained. Assuming that after recalculation, the gas phase volume fraction is reduced by a certain proportion relative to the original result, record this change amount.

[0152] The compensation model assigns weights to three factors: temperature, pressure, and solid phase concentration, which are assumed to be 0.3, 0.4, and 0.3 respectively. The change amounts of the gas phase volume fraction caused by the above three factors (the increase amount caused by temperature, the decrease amount caused by pressure, and the change amount caused by solid phase concentration) are weighted and calculated with the original gas phase volume fraction: the corrected gas phase volume fraction = the original gas phase volume fraction + the change amount affected by temperature×0.3 + the change amount affected by pressure×0.4 + the change amount affected by solid phase concentration×0.3, completing the optimization and adjustment of the calculation result in step 502.

[0153] By comprehensively considering bubble motion parameters, acoustic attenuation characteristics, and environmental factors, the local gas void fraction (VF) of each grid cell can be accurately calculated. Compared to traditional methods that rely solely on total volume changes to estimate VF, this technology can capture subtle gas intrusion events and promptly detect even small amounts of gas entering the drilling fluid, improving the sensitivity and accuracy of VF monitoring. Accurate local VF distribution calculation provides a reliable data foundation for early warning of gas intrusion. By tracking VF changes in each grid cell in real time, the onset and diffusion trend of abnormal VF changes can be quickly identified, shortening the response time from signal acquisition to warning output compared to existing technologies. Once signs of gas intrusion are detected, an early warning can be triggered, buying valuable time for well control measures and effectively reducing the risk of gas intrusion accidents. By considering the impact of environmental factors such as drilling fluid temperature, pressure, and solids concentration on VF and applying compensation corrections, the calculated results are more consistent with actual downhole conditions. Detailed local VF distribution information facilitates in-depth analysis of gas-liquid two-phase flow characteristics within the riser, providing data support for optimizing the drilling fluid circulation system and adjusting drilling parameters. For example, the drilling pump displacement can be adjusted according to the gas content distribution to maintain stable annular pressure, avoid accidents such as well kicks and blowouts caused by gas invasion, and improve drilling operation efficiency and safety.

[0154] In a preferred embodiment of the present invention, step 6, which dynamically tracks the local gas fraction distribution over time, identifies the starting position, diffusion trend, and acceleration of abnormal changes in gas fraction, generates an early warning signal for gas intrusion, and performs real-time pattern matching between the warning signal and a preset threshold to automatically trigger well control commands in a hierarchical manner, may include:

[0155] Step 600 , continuously caching the local gas fraction distribution values of the current and previous two monitoring periods, and constructing a spatiotemporal distribution matrix of gas fraction, where the matrix dimensions include grid coordinates, gas fraction values, and timestamps;

[0156] Step 601: In the spatiotemporal evolution matrix, detect areas where the gas fraction increment of the same grid cell in adjacent periods is greater than or equal to a first threshold, and mark them as abnormal source cells;

[0157] Step 602, taking the abnormal source unit as the center, analyzes the gas content gradient change direction of adjacent grids in the axial or radial direction, and predicts the diffusion vector of the abnormal area in combination with the bubble motion vector set, specifically including:

[0158] Step 6020: extract the three-dimensional grid coordinates of the abnormal source unit, and obtain the local gas content distribution values of the upstream grid unit and the downstream grid unit adjacent to the abnormal source unit in the axial direction, and the circumferential grid unit adjacent to the abnormal source unit in the radial direction;

[0159] Step 6021: Based on the local gas fraction distribution value, calculate the gas fraction difference between the abnormal source unit and each adjacent unit to generate gas fraction gradient change direction data;

[0160] Step 6022: According to the average bubble motion velocity vectors of the grid where the abnormal source unit is located and its adjacent grids in the bubble motion vector set, perform vector superposition correction on the data of the gas holdup gradient change direction to generate an abnormal area diffusion vector.

[0161] Step 603: Calculate the acceleration of the gas holdup change in the abnormal area. When the acceleration is greater than or equal to the second threshold, generate a first-level warning signal; when the abnormal area diffusion vector points to the wellhead direction and the cumulative increase in gas holdup is greater than or equal to the third threshold, generate a second-level warning signal.

[0162] Step 604: Rank the warning signals by level and automatically trigger the corresponding well control instructions, including automatically adjusting the displacement of the drilling pump to compensate for the annulus pressure for the first-level warning, and triggering the pre-closure of the half-closed ram of the blowout preventer group and pressurizing the choke manifold for the second-level warning.

[0163] In the embodiment of the present invention, continuously monitor the local gas holdup distribution values of each grid unit in the riser. Taking each monitoring period (assumed to be 10 seconds) as a time node, cache the gas holdup data of the current period and the previous two periods (a total of three periods, that is, within 30 seconds). Each gas holdup data corresponds to a specific grid unit and includes three-dimensional grid coordinates (axial position, circumferential angle, radial distance), the specific value of the gas holdup, and the timestamp for collecting this data. For example, in the current period, for the grid unit located at an axial position of 20 meters, a circumferential angle of 45°, and a radial distance of 0.3 meters, the gas holdup is 5.5× , and the timestamp is recorded as 10:00:10 on October 1, 2024; integrate the information of all grid units within the three periods and arrange them in the order of grid coordinates, gas holdup values, and timestamps to construct a multi-dimensional gas holdup spatio-temporal distribution matrix. This matrix is like a "data warehouse" that stores the spatio-temporal change information of the gas holdup in the riser within a certain time range.

[0164] Step 601: In the constructed gas holdup spatio-temporal distribution matrix, analyze each grid unit. Compare the gas holdup values of the same grid unit in two adjacent monitoring periods and calculate the gas holdup increment. For example, for a certain grid unit, the gas holdup is 3× in the first period and becomes 5× in the second period, then the gas holdup increment is 5× - 3× = 2× . Compare the calculated gas holdup increment with a pre-set first threshold (assumed to be 1.5× Compare. When the increment of gas holdup is greater than or equal to the first threshold, it indicates that there is a significant change in the gas holdup of this grid cell, and this grid cell is marked as an abnormal source cell. By traversing all grid cells in the matrix, all abnormal source cells that meet the conditions are found, and these cells may be the starting positions of gas invasion.

[0165] Step 6020, for each marked abnormal source cell, first extract its three-dimensional grid coordinates. Taking the abnormal source cell with an axial length of 20 m, a circumferential angle of 45°, and a radial distance of 0.3 m as an example, in the axial direction, determine its upstream grid cell (assuming an axial length of 19 m, a circumferential angle of 45°, and a radial distance of 0.3 m) and downstream grid cell (assuming an axial length of 21 m, a circumferential angle of 45°, and a radial distance of 0.3 m); in the radial direction, obtain its circumferentially adjacent grid cells (such as grid cells with a circumferential angle of 30° and a radial distance of 0.3 m and a circumferential angle of 60° and a radial distance of 0.3 m). Then, obtain the local gas holdup distribution values of these adjacent grid cells at the same moment from the gas holdup spatio-temporal distribution matrix.

[0166] Step 6021, after obtaining the gas holdup distribution values of the abnormal source cell and its adjacent cells, calculate the gas holdup difference between the abnormal source cell and each adjacent cell. For example, the gas holdup of the abnormal source cell is 5× , and the gas holdup of its upstream adjacent cell in the axial direction is 3× , then the gas holdup difference between the two is 5× -3× =2× ; the gas holdup of the downstream adjacent cell in the axial direction is 6× , and the difference is 6× -5× =1× . According to the positive / negative and magnitude of the gas holdup difference, determine the direction of gas holdup gradient change. If the difference from the upstream cell is positive, it indicates that the gas holdup increases upward along the axial direction; if the difference from the downstream cell is negative, it indicates that the gas holdup decreases downstream along the axial direction. Organize this direction information into gas holdup gradient change direction data to describe the change trend of gas holdup in the abnormal area.

[0167] Step 6022: Refer to the bubble motion vector set to obtain the average bubble motion velocity vectors of the grid where the abnormal source unit is located and its adjacent grids. The average bubble motion velocity vector contains the direction and velocity information of the bubble motion. For example, the direction of the average bubble motion velocity vector of the abnormal source unit is upward along the axis, and the velocity is 0.05 m / s. The adjacent units also have their corresponding vectors. Perform vector superposition correction on the gas holdup gradient change direction data and the average bubble motion velocity vector. If the gas holdup gradient change direction is consistent with the bubble motion direction, enhance the trend in this direction; if the directions are inconsistent, make a comprehensive adjustment according to their weights. For example, if the gas holdup gradient shows an increase upward along the axis and the bubble motion direction is also upward along the axis, then the diffusion trend in this direction will be strengthened. In this way, finally generate an abnormal area diffusion vector that can accurately reflect the diffusion direction and velocity of the abnormal area.

[0168] Step 603: Calculate the acceleration by comparing the changes in the gas holdup change amounts in adjacent monitoring periods. For example, the gas holdup change amount from the first period to the second period is 2× , and the gas holdup change amount from the second period to the third period is 3× . Then, the acceleration of the gas holdup change can be obtained by calculating the relationship between the difference in the change amounts (3× - 2× ) and the time interval. Compare the calculated acceleration with the second threshold (assumed to be 0.5× / s²). When the acceleration is greater than or equal to the second threshold, it indicates that the gas holdup change speed is accelerating, and the gas invasion situation may deteriorate rapidly. At this time, generate a first-level warning signal. At the same time, check the direction of the abnormal area diffusion vector and the cumulative increase in gas holdup. When the abnormal area diffusion vector points in the direction of the wellhead and the cumulative increase in gas holdup (i.e., the total increase in gas holdup from the start of the abnormality to the current moment) is greater than or equal to the third threshold (assumed to be 3× ), it means that there is a risk of rapid diffusion of gas invasion towards the wellhead, which may cause serious well control problems. At this time, generate a second-level warning signal.

[0169] Step 604: Rank the generated warning signals. The rank of the secondary warning signal is higher than that of the primary warning signal. According to the rank of the warning signal, the corresponding well control instructions are automatically triggered. When the primary warning signal is received, the displacement of the drilling pump is automatically adjusted. By increasing or decreasing the displacement of the drilling pump, the annulus pressure is compensated to prevent the further development of gas invasion. For example, if it is detected that gas invasion causes a decrease in the annulus pressure, the system will automatically increase the displacement of the drilling pump to increase the annulus pressure and maintain the downhole pressure balance. When the secondary warning signal is received, the situation is more urgent. Immediately trigger the pre - closing of the half - closed ram of the blowout preventer group and pressurize the choke manifold at the same time. The pre - closing of the half - closed ram can, to a certain extent, prevent the gas from continuing to diffuse upward, and pressurizing the choke manifold helps control the wellhead pressure and prevent blowout accidents, ensuring the safety of drilling operations to the greatest extent.

[0170] By constructing a spatio - temporal distribution matrix of gas holdup and dynamically tracking the local gas holdup through time series, it is possible to promptly capture the minute abnormal changes in gas holdup and accurately identify the starting position of gas invasion. Compared with traditional monitoring methods, it can issue a warning when gas invasion just occurs and the change in gas holdup is not obvious, greatly advancing the warning time and winning precious time for taking measures to control gas invasion. Using the gradient change of gas holdup and the bubble motion vector to predict the diffusion vector of the abnormal area can not only judge the diffusion direction of gas invasion, but also estimate the diffusion speed and acceleration. This enables drilling workers to understand the development trend of gas invasion in advance, make preparations in advance, and avoid serious accidents such as well kick and blowout caused by the diffusion of gas invasion, ensuring the safety of drilling equipment and personnel.

[0171] Automatically generating warning signals of different levels according to the severity of gas invasion and correspondingly triggering different well control instructions realizes the automation and precision of well control operations. The primary warning compensates the pressure by adjusting the displacement of the drilling pump, and the secondary warning directly starts key equipment such as the blowout preventer and the choke manifold. This hierarchical response mechanism can quickly take the most appropriate measures according to the actual situation, improve the well control efficiency, reduce the probability of accidents, and at the same time reduce the mistakes that may be brought by manual intervention, enhancing the safety and reliability of drilling operations. The detailed gas invasion warning information and well control instruction triggering mechanism provide strong support for drilling operation decision - making. Workers can understand the downhole gas invasion situation according to the warning signal, reasonably adjust drilling parameters, and optimize the operation process. For example, when the gas invasion risk is high, adjust the performance of the drilling fluid in advance and arrange personnel to make emergency preparations to ensure the efficient progress of drilling operations under the premise of safety, reducing operation costs and risks.

[0172] As Figure 2 shown, the embodiment of the present invention also provides an underwater ultrasonic Doppler real - time analysis device based on a digital transducer, including:

[0173] A signal acquisition unit, which is used to transmit ultrasonic waves through a digital transducer array on the outer wall of the riser and receive echo signals, and output an original digital signal containing phase and amplitude information;

[0174] A frequency shift analysis unit, which is used to extract Doppler frequency shift characteristics from the original digital signal and output the real-time movement speed and direction of the bubbles;

[0175] A grid mapping unit, which is used to establish a three-dimensional dynamic space grid and discretely map the bubble movement parameters, and output a set of bubble movement vectors for each grid unit;

[0176] A vector correction unit, which is used to perform gradient analysis and dynamic correction on the set of bubble movement vectors of adjacent grid units, and output optimized bubble movement parameters;

[0177] A gas holdup calculation unit, which is used to calculate the local gas holdup distribution in each grid unit in real time based on the optimized bubble movement parameters and combined with the analysis of the acoustic attenuation characteristics on the ultrasonic wave propagation path;

[0178] An early warning execution unit, which is used to dynamically track the local gas holdup distribution and trigger well control instructions in a graded manner.

[0179] It should be noted that this device corresponds to the above method. All implementation manners in the above method embodiments are applicable to this embodiment and can also achieve the same technical effects.

[0180] The above is the preferred implementation manner of the present invention. It should be pointed out that for those of ordinary skill in the art, without departing from the principle described in the present invention, several improvements and refinements can be made, and these improvements and refinements should also be regarded as the protection scope of the present invention.

Claims

1. An underwater ultrasonic Doppler real-time analysis method based on digital transducers, characterized in that, The method includes the following steps: Step 1: Transmit multiple beams of directional ultrasonic waves into the drilling fluid through a digital transducer array on the outer wall of the riser, receive the echo signals reflected by the gas-liquid two-phase flow, and perform real-time digital conversion on the echo signals to generate original digital signals containing phase and amplitude information; Step 2: Extract the Doppler frequency shift characteristics of the original digital signals, and determine the real-time movement speed and direction of the bubbles in the drilling fluid by calculating the frequency deviation between the reflected ultrasonic waves and the transmitted ultrasonic waves; Step 3: Establish a three-dimensional dynamic spatial grid based on the flow field characteristics inside the riser, and perform spatial discretization mapping of the real-time movement speed and direction of the bubbles according to the grid positions to generate a set of bubble movement vectors for each grid cell; Step 4: Analyze the gradient changes of the set of bubble movement vectors in adjacent grid cells, generate grid characteristic difference values according to the sudden change of vector direction and rate difference, and dynamically correct the set of bubble movement vectors based on the grid characteristic difference values, and output the optimized bubble movement parameters; Step 5: Based on the optimized bubble movement parameters, combined with the analysis of the acoustic attenuation characteristics on the ultrasonic wave propagation path, calculate the local gas holdup distribution in each grid cell in real time; Step 6: Perform time-series dynamic tracking on the local gas holdup distribution, identify the starting position, diffusion trend and acceleration of abnormal gas holdup changes, generate early warning signals for gas invasion, and perform real-time pattern matching of the warning signals with preset thresholds to automatically trigger well control instructions at different levels.

2. The real-time underwater ultrasonic Doppler analysis method based on a digital transducer according to claim 1, wherein Transmitting multiple beams of directional ultrasonic waves into the drilling fluid through a digital transducer array on the outer wall of the riser, receiving the echo signals reflected by the gas-liquid two-phase flow, and performing real-time digital conversion on the echo signals to generate original digital signals containing phase and amplitude information includes: Capturing ultrasonic echo analog signals at different axial depth positions through a transducer array arranged annularly on the outer wall of the riser; Detecting the spectral characteristics of the drilling fluid background noise, removing the components matching the noise spectrum from the ultrasonic echo analog signals to generate primary purification signals; According to the change of the acoustic characteristics of the drilling fluid in the deep-water high-pressure environment, adjusting the gain amplitude of the primary purification signals in real time to generate secondary optimized signals with balanced intensity; Discretely sampling the secondary optimized signals at preset time intervals, recording the amplitude values, phase angles and accurate timestamps of each sampling point to form a digital signal sequence; Performing phase alignment and amplitude superposition on multi-channel digital signals with overlapping time windows and spatial coordinates to generate original digital signals with three-dimensional position tags.

3. The underwater ultrasonic Doppler real-time analysis method based on a digital transducer according to claim 2, wherein, Extracting the Doppler frequency shift characteristics of the original digital signals, and determining the real-time movement speed and direction of the bubbles in the drilling fluid by calculating the frequency deviation between the reflected ultrasonic waves and the transmitted ultrasonic waves includes: Dividing the original digital signals with three-dimensional position tags into continuous time windows at a fixed duration, parsing the ultrasonic frequency offset and phase change trajectory in each window, and outputting a time-frequency feature matrix; Based on the time-frequency feature matrix, identifying signal regions with continuously similar frequency shift amplitudes and consistent phase evolution directions, and dividing them into independent bubble acoustic fingerprint clusters; For each cluster of bubble acoustic fingerprints, calculate the difference in frequency shift amount between adjacent time windows along the ultrasonic emission direction to generate the axial motion rate component of the bubble; calculate the tangential offset angle component of the bubble group based on the phase difference between adjacent receiving units of the bubble acoustic fingerprint cluster.

4. The underwater ultrasonic Doppler real-time analysis method based on a digital transducer according to claim 3, characterized in that Based on the time-frequency feature matrix, identify signal regions with continuously similar frequency shift amplitudes and consistent phase evolution directions, and segment them into independent bubble acoustic fingerprint clusters, including: In the time-frequency feature matrix, extract the frequency shift amplitude sequence of the same spatial coordinate point within three consecutive time windows, calculate the absolute difference in frequency shift amount between adjacent windows, and if the absolute differences are all less than the dynamic fluctuation threshold, mark this coordinate point as a frequency shift stable point; For the cluster of frequency shift stable points, analyze the phase change trajectory, calculate the phase difference symbol combination between adjacent receiving units within the same time window, and when the phase difference symbol combination meets the preset direction consistency condition, determine that the phase evolution direction is consistent; Cluster the cluster of frequency shift stable points that meet the consistent phase evolution direction into independent bubble acoustic fingerprint clusters according to the principle of spatial continuity.

5. The underwater ultrasonic Doppler real-time analysis method based on a digital transducer according to claim 4, characterized in that, Based on the flow field characteristics in the riser, establish a three-dimensional dynamic space grid, spatially discretize and map the real-time motion speed and direction of the bubbles according to the grid positions to generate the bubble motion vector set of each grid cell, including: Fuse the axial motion rate component and the tangential offset angle component of the bubble to generate a three-dimensional bubble motion vector with spatial coordinate labels; According to the inner diameter size of the riser and the characteristics of the gas-liquid two-phase flow, divide the monitoring area into annular grid cells that are stratified axially and sectorized circumferentially, and distribute the three-dimensional bubble motion vectors to the corresponding annular grid cells according to the spatial coordinate labels; Perform timestamp verification on the vector data that persists for a preset duration in each grid cell to generate the current valid bubble motion vector set.

6. The underwater ultrasonic Doppler real-time analysis method based on a digital transducer according to claim 5, characterized in that, Perform gradient change analysis on the bubble motion vector sets of adjacent grid cells, generate the grid characteristic difference value based on the sudden change of the vector direction and the rate difference, and dynamically correct the bubble motion vector set based on the grid characteristic difference value, and output the optimized bubble motion parameters, including: For the current valid bubble motion vector set, calculate the maximum deviation value and the rate change rate of the motion direction angle between each grid cell and its adjacent cells; Generate the grid characteristic difference value based on the maximum deviation value and the rate change rate; Compare the grid characteristic difference value with the preset threshold. When the grid characteristic difference value is greater than the preset threshold, generate a strong correction factor; when the grid characteristic difference value is less than the preset threshold, generate a weak correction factor; Based on the strong correction factor and the weak correction factor, perform neighborhood smoothing on the motion vector of the current grid cell, and output the optimized bubble motion parameters.

7. The underwater ultrasonic Doppler real-time analysis method based on a digital transducer according to claim 6, characterized in that Based on the optimized bubble motion parameters, combined with the analysis of the acoustic attenuation characteristics on the ultrasonic propagation path, calculate the local gas holdup distribution in each grid cell in real time, including: According to the optimized bubble motion parameters, extract the average motion rate of the bubbles in each grid cell and the number of bubbles per unit volume, and calculate the sound intensity attenuation rate before and after the ultrasonic wave penetrates the current grid cell; Calculate the additional acoustic attenuation coefficient caused by the bubble group according to the sound intensity attenuation rate and the number of bubbles. Determine the gas volume fraction through the correlation between the additional sound attenuation coefficient and the equivalent radius of the bubbles; Collect the temperature, pressure, and solid concentration of the drilling fluid in real time, compensate for the gas volume fraction, and output the local gas holdup distribution values with three-dimensional grid coordinates.

8. The underwater ultrasonic Doppler real-time analysis method based on a digital transducer according to claim 7, characterized in that Conduct time-series dynamic tracking on the local gas holdup distribution, identify the starting position, diffusion trend, and acceleration of abnormal changes in the gas holdup, generate early warning signals for gas invasion, and perform real-time pattern matching of the warning signals with preset thresholds to automatically trigger well control commands at different levels, including: Continuously cache the local gas holdup distribution values of the current and the previous two monitoring cycles, and construct a spatio-temporal distribution matrix of the gas holdup. The matrix dimensions include grid coordinates, gas holdup values, and timestamps; In the spatio-temporal evolution matrix, detect the areas where the gas holdup increment of the same grid cell in adjacent cycles is greater than or equal to the first threshold, and mark them as abnormal source cells; Taking the abnormal source cells as the center, analyze the direction of the gas holdup gradient change in the axial or radial direction of adjacent grids, and combine the bubble motion vector set to predict the diffusion vector of the abnormal area; Calculate the acceleration of the gas holdup change in the abnormal area. When the acceleration is greater than or equal to the second threshold, generate a first-level warning signal; when the diffusion vector of the abnormal area points towards the wellhead direction and the cumulative increase in the gas holdup is greater than or equal to the third threshold, generate a second-level warning signal; Rank the warning signals by level and automatically trigger corresponding well control commands, including automatically adjusting the displacement of the drilling pump to compensate for the annulus pressure for the first-level warning, and triggering the pre-closure of the half-closed ram of the blowout preventer group and pressurizing the choke manifold for the second-level warning.

9. The underwater ultrasonic Doppler real-time analysis method based on a digital transducer according to claim 8, characterized in that Taking the abnormal source cells as the center, analyze the direction of the gas holdup gradient change in the axial or radial direction of adjacent grids, and combine the bubble motion vector set to predict the diffusion vector of the abnormal area, including: Extract the three-dimensional grid coordinates of the abnormal source cells, and obtain the local gas holdup distribution values of the upstream grid cell and the downstream grid cell adjacent axially to the abnormal source cells, and the circumferential grid cells adjacent radially; Based on the local gas holdup distribution values, calculate the gas holdup difference between the abnormal source cells and each adjacent cell, and generate data on the direction of the gas holdup gradient change; According to the average bubble motion velocity vectors of the grids where the abnormal source cells are located and the adjacent grids in the bubble motion vector set, perform vector superposition correction on the data on the direction of the gas holdup gradient change to generate the diffusion vector of the abnormal area.

10. An underwater ultrasonic Doppler real-time analysis device based on a digital transducer, which implements the method according to any one of claims 1 to 9, characterized in that, Including: A signal acquisition unit for transmitting ultrasonic waves through a digital transducer array on the outer wall of the riser and receiving echo signals, and outputting the original digital signals containing phase and amplitude information; A frequency shift analysis unit for extracting the Doppler frequency shift characteristics of the original digital signals and outputting the real-time motion velocity and direction of the bubbles; A grid mapping unit for establishing a three-dimensional dynamic space grid and discretely mapping the bubble motion parameters, and outputting the bubble motion vector set of each grid cell; A vector correction unit for performing gradient analysis and dynamic correction on the bubble motion vector sets of adjacent grid cells, and outputting the optimized bubble motion parameters; A gas holdup calculation unit for calculating the local gas holdup distribution in each grid cell in real time based on the optimized bubble motion parameters and combining the acoustic attenuation characteristics analysis on the ultrasonic wave propagation path; An early warning execution unit is used to dynamically track the local gas holdup distribution and trigger well control instructions in a graded manner.

Images