On-line correction and metering method for ultrasonic natural gas flowmeter

Abstract

The present application relates to natural gas flow metering technical field, specifically disclose an ultrasonic natural gas flow meter on-line correction and metering method, collecting each sound channel propagation time and solving flow velocity and measured sound velocity, and obtaining component and temperature and pressure parameters; the measured sound velocity of each sound channel is arranged along the radial direction, and the low-frequency gradient component and the high-frequency random component are separated out through the variational mode decomposition; the radial monotonicity index is calculated according to the low-frequency gradient component, which is matched with the chromatographic component change time stamp in time window, and the disturbance source identification distinguishing component disturbance and hardware drift is output; when the component disturbance is kept, the sound path length and the timing reference are unchanged, the sound velocity deviation is converted into the fine adjustment amount of compression factor, and when the hardware drift is corrected, the sound path length and the timing reference are kept unchanged, and the compression factor is kept unchanged; the cross-sectional average flow velocity is reconstructed according to the corrected parameters, and the standard reference volume flow is converted and output; the present application can distinguish different disturbance sources and correct respectively, and avoid the measurement deviation caused by false correction.

Description

Technical Field

[0001] This invention relates to the field of natural gas flow metering technology, specifically to an online correction and metering method for ultrasonic natural gas flow meters. Background Technology

[0002] Ultrasonic flow meters, due to their advantages such as having no moving parts, no pressure loss, and a wide range, have been widely used in the field of natural gas trade and transfer metering. Their basic principle is the time-of-flight method: by measuring the time difference between the downstream and upstream propagation of ultrasonic pulses in the medium, and combining this with a fixed sound path length, the average velocity along the sound channel is calculated. Then, through multi-channel arrangement and numerical integration methods, the average velocity of the pipe cross-section is reconstructed, and finally converted into the volumetric flow rate or energy flow rate under standard reference conditions.

[0003] In existing online correction methods for ultrasonic natural gas flow meters, the state component introduced by the non-uniform mixing of medium components and the motion component introduced by the thermal expansion and contraction of the transducer and the drift of the timing reference in the measured sound velocity deviation are not distinguished. This leads to the misjudgment of the radial sound velocity gradient change caused by component stratification as hardware parameter drift, which in turn leads to the reverse correction of the sound path length and timing reference. As a result, the deviation caused by component fluctuation is solidified into the hardware measurement reference, forming an irreversible measurement offset. Summary of the Invention

[0004] The purpose of this invention is to provide an online correction and metering method for ultrasonic natural gas flow meters to solve the problems mentioned above.

[0005] The objective of this invention can be achieved through the following technical solutions: The online calibration and metering method for ultrasonic natural gas flow meters includes the following steps: S1, collect the propagation time of ultrasonic waves in both directions of flow for each channel, calculate the line average flow velocity and the line average measured sound velocity for each channel, and simultaneously obtain the mole fraction of natural gas components and pipeline temperature and pressure parameters; S2, the measured sound velocities of each channel are arranged radially along the pipe, and variational mode decomposition is performed on them to separate the low-frequency large-scale gradient component that characterizes the component stratification effect and the high-frequency narrowband random component that characterizes the thermal expansion and contraction of the transducer and the drift of the timing reference. S3 calculates the radial monotonicity index based on the low-frequency large-scale gradient components, and matches the monotonicity index with the time stamp of the upstream chromatographic component change to output a disturbance source identifier to distinguish between state component disturbance and motion hardware drift. S4, when the disturbance source indicator indicates component disturbance, keep the current values ​​of sound path length and timing reference unchanged, and only convert the deviation between the measured sound velocity and the theoretical sound velocity into the compression factor fine adjustment amount according to the preset mapping relationship; when the indicator indicates hardware drift, activate the correction of sound path length and timing reference, keep the compression factor unchanged, and output the corrected sound path length, timing reference and current effective compression factor. S5, based on the corrected sound path length, timing reference and current effective compression factor, reconstructs the cross-sectional average flow velocity by combining the average flow velocity of each channel line, and converts the working condition volumetric flow rate into the standard volumetric flow rate under standard reference conditions, as the output of the measurement result.

[0006] As a further aspect of the present invention: S2 specifically includes: The measured sound velocities of each channel are used as a spatial sequence that is discretely distributed along the radial direction of the pipe. The Hilbert transform is then applied to the spatial sequence to obtain the analytical signals of each modal component. The center frequency and bandwidth of each modal component are updated iteratively, and in each iteration, the upper limit of the bandwidth of the low-frequency large-scale gradient component is constrained to be less than the first preset value and its radial second-order difference norm is constrained to be less than the second preset value. At the same time, the lower limit of the bandwidth of the high-frequency narrowband random component is constrained to be greater than the third preset value. The iteration terminates when the sum of the center frequency changes of all modal components in two adjacent iterations is less than the fourth preset value, and the separated low-frequency large-scale gradient component and high-frequency narrowband random component are output.

[0007] As a further aspect of the present invention: the analytic signals of each modal component obtained after performing a Hilbert transform on the spatial sequence specifically include: Using the measured sound velocity of each channel as the initial real part of the analytical signal, the measured sound velocity distribution sequence along the radial direction is extended to generate a radially extended sequence. The Hilbert transform is applied to the radially extended sequence to generate an imaginary part sequence orthogonal to the real part. The real part sequence and the imaginary part sequence together constitute the analytical signal of each modal component. Extract the amplitude envelope corresponding to the original radial position in the analytical signal, perform differential operation between adjacent channels on the amplitude envelope, and output the amplitude change between each adjacent channel; The amplitude changes between each adjacent channel are compared with the preset radial continuity tolerance. When all amplitude changes are less than the radial continuity tolerance, the current analytical signal is output to the iterative update stage to participate in subsequent iterative calculations. Otherwise, radial smoothing is applied to the imaginary part sequence of the current analytical signal until the smoothed amplitude changes meet the radial continuity tolerance before outputting.

[0008] As a further aspect of the present invention: S3 specifically includes: For the low-frequency large-scale gradient component, the sound velocity difference between adjacent channels is calculated sequentially along the radial direction. The sign value of each sound velocity difference is accumulated, and the ratio of the absolute value of the accumulated result to the total number of channels is used as the monotonicity index output. Extract the moment when the mole fraction of a component jumps from the upstream chromatographic data as the component change timestamp, calculate the difference between the current time and the component change timestamp, and mark the time window as successfully matched when the absolute value of the difference is less than the preset time window value. When the monotonicity exponent is greater than the preset monotonicity threshold and the time window is successfully matched, the disturbance source identifier is set to the component disturbance state; otherwise, it is set to the hardware drift state, and the disturbance source identifier is output.

[0009] As a further aspect of the present invention: S4 specifically includes: The type of disturbance is determined based on the source identifier. When the identifier is a component disturbance, the measured sound velocities of each channel are weighted and averaged radially, and the difference is calculated with the theoretical sound velocity to obtain the sound velocity deviation. The corresponding compression factor fine-tuning amount is determined by piecewise linear interpolation based on the absolute value of the sound velocity deviation. The compression factor fine-tuning amount is then superimposed on the theoretical compression factor to output the current effective compression factor, while keeping the sound path length and timing reference at their current values. When hardware drift is identified, the average deviation between the measured sound velocity and the theoretical sound velocity of each channel is used as the correction basis. The average deviation is proportionally calculated and used as the path length correction amount and the timing reference correction amount, respectively. These are then superimposed on the current path length and the current timing reference, and the corrected path length and the corrected timing reference are output. At the same time, the theoretical compression factor is kept as the current effective compression factor. The current effective compression factor, the corrected path length, and the corrected timing reference are all used as input parameters for S5.

[0010] As a further aspect of the present invention: the determination of the corresponding compression factor fine-tuning amount after piecewise linear interpolation specifically includes: Multiple continuous numerical intervals between the absolute value of the sound velocity deviation and the compression factor fine-tuning amount are pre-constructed, and corresponding interval start fine-tuning value and interval end fine-tuning value are configured for each numerical interval. The absolute value of the current sound speed deviation is compared with the boundary value of each numerical interval to locate the target numerical interval into which the absolute value falls, and the relative position percentage of the absolute value within the target numerical interval is calculated. Based on the relative position ratio, interpolate linearly between the starting fine-tuning value and the ending fine-tuning value of the target value interval, output the compression factor fine-tuning amount corresponding to the absolute value, and superimpose the compression factor fine-tuning amount to the theoretical compression factor to obtain the current effective compression factor.

[0011] As a further aspect of the present invention: the conversion of the operating condition volumetric flow rate into the standard volumetric flow rate under standard reference conditions, and outputting it as the measurement result, specifically includes: Based on the corrected path length, the path deviation of the line average velocity corresponding to each channel is corrected to obtain the corrected velocity of each channel. Then, the timing deviation is inverted and corrected based on the time difference between the forward and reverse flow corresponding to the corrected velocity of each channel using the corrected timing reference to obtain the final velocity of each channel. Each channel is assigned a preset radial position weighting coefficient, and the sum of the products of the final flow velocity of each channel and its radial position weighting coefficient is used as the cross-sectional average flow velocity. Using the current effective compressibility factor as the compressibility factor value under standard reference conditions, the operating condition volumetric flow rate is converted to the standard reference volumetric flow rate at the standard reference pressure and standard reference temperature according to the pressure ratio and temperature ratio, and the standard reference volumetric flow rate is output as the measurement result.

[0012] As a further aspect of the present invention: obtaining the final flow velocity of each channel specifically includes: Obtain the deviation ratio between the corrected timing reference and the standard timing reference. Perform an inverse proportional operation on the original forward and reverse flow time differences corresponding to each channel with the deviation ratio to obtain the corrected forward flow time difference and the corrected reverse flow time difference. The difference between the corrected downstream time difference and the corrected upstream time difference is taken as the effective time difference. The effective time difference is multiplied by the corrected sound path length and then divided by the current measured sound speed to obtain the final time difference compensation amount of the sound channel. The final time difference compensation is superimposed on the corrected downstream time difference and the corrected upstream time difference. The superimposed downstream time difference and upstream time difference are used as the basis for calculating the final flow rate of each channel, and the final flow rate of each channel is output.

[0013] The beneficial effects of this invention are: (1) By performing variational mode decomposition on the measured sound velocity of each channel and matching the monotonicity index with the chromatographic timestamp to output the disturbance source identifier, it is possible to clearly distinguish the sound velocity change caused by component stratification and the sound velocity change caused by transducer hardware drift. This allows for the separate execution of compression factor fine-tuning or hardware parameter correction for different disturbance sources, avoiding the misjudgment of component fluctuations as hardware drift and the reverse correction of sound path length and timing reference, thereby preventing irreversible metrological offset caused by miscorrection.

[0014] (2) When the component disturbance is determined, the sound path length and timing reference remain unchanged, and the sound velocity deviation is mapped to the compression factor fine adjustment amount, so that the state quantity deviation caused by the non-uniform mixing of components is limited to the compression factor correction channel and does not contaminate the hardware measurement reference channel; when the hardware drift is determined, only the sound path length and timing reference are corrected, without changing the compression factor, thus realizing the physical isolation of the correction paths of the two types of physical quantities, ensuring the integrity of metrological traceability and the traceability of correction values ​​during the online correction process. Attached Figure Description

[0015] The invention will now be further described with reference to the accompanying drawings.

[0016] Figure 1 This is a flowchart of the method of the present invention; Figure 2 This is a flowchart of the process of obtaining the analytical signals of each modal component in this invention; Figure 3 This is a flowchart illustrating the process of determining the corresponding compression factor fine-tuning amount after piecewise linear interpolation in this invention. Detailed Implementation

[0017] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0018] Please see Figure 1 As shown, this invention provides an online correction and metering method for ultrasonic natural gas flow meters, comprising the following steps: S1, collect the propagation time of ultrasonic waves in both directions of flow for each channel, calculate the line average flow velocity and the line average measured sound velocity for each channel, and simultaneously obtain the mole fraction of natural gas components and pipeline temperature and pressure parameters; S2, the measured sound velocities of each channel are arranged radially along the pipe, and variational mode decomposition is performed on them to separate the low-frequency large-scale gradient component that characterizes the component stratification effect and the high-frequency narrowband random component that characterizes the thermal expansion and contraction of the transducer and the drift of the timing reference. S3 calculates the radial monotonicity index based on the low-frequency large-scale gradient components, and matches the monotonicity index with the time stamp of the upstream chromatographic component change to output a disturbance source identifier to distinguish between state component disturbance and motion hardware drift. S4, when the disturbance source indicator indicates component disturbance, keep the current values ​​of sound path length and timing reference unchanged, and only convert the deviation between the measured sound velocity and the theoretical sound velocity into the compression factor fine adjustment amount according to the preset mapping relationship; when the indicator indicates hardware drift, activate the correction of sound path length and timing reference, keep the compression factor unchanged, and output the corrected sound path length, timing reference and current effective compression factor. S5, based on the corrected sound path length, timing reference and current effective compression factor, reconstructs the cross-sectional average flow velocity by combining the average flow velocity of each channel line, and converts the working condition volumetric flow rate into the standard volumetric flow rate under standard reference conditions, as the output of the measurement result.

[0019] In S1, the propagation time of ultrasonic waves in both upstream and downstream directions for each channel is collected, and the line-average flow velocity and the line-average measured sound velocity for each channel are calculated. Simultaneously, the mole fraction of natural gas components and pipeline temperature and pressure parameters are obtained, specifically including: Multiple pairs of ultrasonic transducers are arranged at different chord heights along the pipe cross-section. Each pair of transducers is installed on opposite side walls of the pipe, forming an ultrasonic channel. Four channels are arranged along the pipe cross-section, located at different radial heights. During each measurement cycle, the two transducers on each channel alternately transmit and receive ultrasonic pulse signals. One transducer transmits ultrasonic waves along the airflow direction, while the other receives them against the airflow direction, obtaining the downstream propagation time. This process is then repeated in reverse to obtain the upstream propagation time. The downstream and upstream propagation times are directly captured and output by a high-precision time-to-digital converter circuit with a timing resolution better than one nanosecond.

[0020] After obtaining the downstream and upstream propagation times for each channel, the downstream and upstream propagation times are summed, and the sum is divided by twice the path length to obtain the line-average measured sound velocity for that channel. The path length refers to the fixed physical distance between the two transducers in that channel, and its value is determined during the flow meter's factory calibration and stored in a non-volatile memory unit.

[0021] The difference between the downstream propagation time and the upstream propagation time is calculated, and then multiplied by twice the sound path length. This difference is then divided by the square of the sum of the downstream and upstream propagation times to obtain the line-average velocity of that channel. Using this method, the line-average measured sound velocity and the line-average velocity of each channel are calculated individually.

[0022] Regarding the acquisition of natural gas component data, an online gas chromatograph is installed at the upstream straight pipe section of the ultrasonic flow meter. The gas chromatograph continuously extracts natural gas samples from the pipeline at a preset sampling cycle. After separation by the chromatographic column, the detector outputs the mole fraction of each component (including methane, ethane, propane, n-butane, isobutane, n-pentane, isopentane, hexane, nitrogen, and carbon dioxide).

[0023] In this embodiment, the sampling cycle of the gas chromatograph is 30 seconds, and after each sampling is completed, the mole fraction of each component is transmitted to the flow calculation unit via digital communication.

[0024] For obtaining the static pressure of the pipeline, a pressure tap is opened in the ultrasonic flow meter body or an adjacent straight pipe section, and a pressure transmitter is installed. This pressure transmitter continuously senses the absolute pressure value of the medium in the pipeline with a millisecond-level response speed and converts the sensed pressure value into a standard current signal output. For obtaining the medium temperature, a platinum resistance temperature transmitter is installed at the same location. This platinum resistance temperature transmitter extends into the pipeline and directly contacts the medium to sense the current temperature value of the medium and converts the temperature value into a standard resistance signal output. The signals from the pressure transmitter and the platinum resistance temperature transmitter are synchronously transmitted to the flow calculation unit after analog-to-digital conversion.

[0025] The propagation time in both directions of flow for each channel, the average measured sound velocity along each channel line, the average flow velocity along each channel line, the mole fraction of natural gas components, the static pressure in the pipeline, and the temperature of the medium are used as the initial input parameters for subsequent steps.

[0026] Among them, the measured linear average sound velocity is arranged radially along the pipe according to the channel number, forming a spatial sequence distributed radially.

[0027] Please see Figure 2 As shown, in S2, the measured sound velocities of each channel are arranged radially along the pipe, and variational mode decomposition is performed on them to separate the low-frequency large-scale gradient components characterizing the component stratification effect and the high-frequency narrowband random components characterizing the thermal expansion and contraction of the transducer and the drift of the timing reference. Specifically, these include: After obtaining the average measured sound velocity of each channel, the average measured sound velocities of each channel are arranged into a spatial sequence distributed radially, according to the radial position of each channel in the pipe cross-section from low to high. The length of this spatial sequence is equal to the total number of channels, and each element in the sequence corresponds to the measured sound velocity value at a radial position.

[0028] Variational mode decomposition is performed on the radial spatial sequence composed of measured sound velocities in each channel.

[0029] The iteration process of variational mode decomposition is as follows: Using the measured sound velocity of each channel as the initial real part of the analytical signal, boundary extensions are performed at both ends of the measured sound velocity distribution sequence along the radial direction. The left extension extends outward by two data points with the difference between the sound velocity values ​​of two adjacent channels at the beginning of the sequence as the extension step size, and the right extension extends outward by two data points with the difference between the sound velocity values ​​of two adjacent channels at the end of the sequence as the extension step size, thus generating a radially extended sequence.

[0030] Applying a Hilbert transform to the radially extended sequence generates an imaginary sequence orthogonal to the real part. The real and imaginary sequences together constitute the analytic signals of each modal component.

[0031] The magnitude of the analytic signal is calculated to obtain the amplitude envelope of the analytic signal, and the interval segment corresponding to the original radial position is extracted from the amplitude envelope.

[0032] The truncated amplitude envelope is subjected to a difference operation between adjacent channels. Specifically, the absolute value of the difference between the amplitude envelope values ​​of two adjacent channels is calculated sequentially, and the amplitude change between each adjacent channel is output. The amplitude change between each adjacent channel is then compared with a preset radial continuity tolerance, which is set to 1% of the average measured sound velocity of each channel.

[0033] When all amplitude changes are less than the radial continuity tolerance, the current analytical signal is output to the iterative update stage; otherwise, radial smoothing is applied to the imaginary part sequence of the current analytical signal. This smoothing process uses a three-point moving average method, that is, the imaginary part value of each data point is replaced by the average of its own and the imaginary part values ​​of the two adjacent data points. The amplitude envelope and difference operation are repeatedly calculated until the smoothed amplitude changes meet the condition of being less than the radial continuity tolerance and then output.

[0034] In the iterative update process, the initial values ​​of the center frequencies of each modal component are first set. The initial value of the center frequency of the low-frequency large-scale gradient component is set to the average value of the absolute values ​​of the first-order difference after the measured sound velocities of all channels are arranged radially. The initial value of the center frequency of the high-frequency narrowband random component is set to twice the above average value.

[0035] Set the initial bandwidth values ​​for each modal component, where the initial bandwidth value for the low-frequency component is set to 5% of the average measured sound velocity of each channel, and the initial bandwidth value for the high-frequency component is set to 1% of the average measured sound velocity of each channel.

[0036] In each iteration, the analytical signal is projected to the vicinity of the center frequency corresponding to the low-frequency large-scale gradient component and the vicinity of the center frequency corresponding to the high-frequency narrowband random component, and each modal component and its corresponding center frequency and bandwidth are updated.

[0037] Constraints are applied during the update process: the upper limit of the bandwidth of the low-frequency large-scale gradient component is less than a first preset value, which is 3% of the average measured sound velocity of each channel; and the second-order difference norm of the low-frequency large-scale gradient component after radial arrangement is less than a second preset value. The second-order difference norm is calculated as follows: calculate the second-order difference value for each of the three adjacent data points in the low-frequency large-scale gradient component, sum the absolute values ​​of each second-order difference value and divide by the total number of data points. The second preset value is 2% of the average measured sound velocity of each channel. The lower limit of the bandwidth of the high-frequency narrowband random component is greater than a third preset value, which is 0.5% of the average measured sound velocity of each channel.

[0038] After each iteration, the change in center frequency of each modal component between the current iteration and the previous iteration is calculated. The changes in center frequency of all modal components are summed to obtain the sum of center frequency changes. The iteration terminates when the sum of center frequency changes is less than a fourth preset value, which is one ten-thousandth of the average measured sound velocity of each channel.

[0039] After the iteration terminates, the low-frequency large-scale gradient component and the high-frequency narrowband random component are output after separation. The low-frequency large-scale gradient component is used as the output to characterize the component stratification effect, and the high-frequency narrowband random component is used as the output to characterize the thermal expansion and contraction of the transducer and the drift of the timing reference.

[0040] In S3, the radial monotonicity index is calculated based on the low-frequency large-scale gradient components, and the monotonicity index is matched with the time window of the upstream chromatographic component change timestamps. This outputs a disturbance source identifier to distinguish between state component disturbances and motion hardware drift, specifically including: After obtaining the low-frequency large-scale gradient component, the sound velocity values ​​of each channel contained in the component are arranged in order from low to high according to the radial position of the channel.

[0041] The sound velocity difference between adjacent channels is calculated sequentially along the radial direction. That is, the sound velocity value at the next position is subtracted from the sound velocity value at the previous position to obtain the sound velocity difference between adjacent channels.

[0042] For each sound speed difference, determine its sign: when the sound speed difference is positive, record its sign value as positive 1; when the sound speed difference is negative, record its sign value as negative 1; when the sound speed difference is 0, record its sign value as numerical 0.

[0043] The sign values ​​corresponding to the sound velocity differences between all adjacent channels are summed to obtain the sign sum result. The absolute value of this sign sum result is divided by the total number of channels, and the resulting ratio is output as the monotonicity index.

[0044] The monotonicity index ranges between 0 and 1. When all sound velocity differences have the same sign, the monotonicity index is 1, indicating that the low-frequency large-scale gradient component exhibits a strictly monotonic distribution along the radial direction.

[0045] In terms of time window matching, records of the change of mole fraction of each component over time are extracted from historical data output by the upstream online gas chromatograph.

[0046] For the component mole fraction data at each sampling time, compare the difference between the current sampling time and the previous sampling time for each component mole fraction. When the absolute value of the difference in mole fraction of any key component (including methane, nitrogen, and carbon dioxide) is greater than the preset jump threshold, mark the current sampling time as the component change timestamp.

[0047] The preset threshold value for the jump is 0.005, meaning that a change in mole fraction exceeding 0.5% is considered a jump.

[0048] If the mole fractions of multiple components change simultaneously at the same time, the time when the change threshold is reached first is taken as the timestamp of the component change at that time.

[0049] Obtain the current measurement time and calculate the difference between it and the most recently obtained component change timestamp to obtain the time interval between the two. Compare the absolute value of this time interval with a preset time window value of three seconds. If the absolute value of the time interval is less than three seconds, the time window is marked as a successful match.

[0050] Regarding the output disturbance source identifier, the matching status of the monotonicity index and the time window is also determined. When the monotonicity index is greater than the preset monotonicity threshold and the time window is successfully matched, the disturbance source identifier is set to the component disturbance state; otherwise, the disturbance source identifier is set to the hardware drift state.

[0051] The preset monotonic threshold is 0.8. The disturbance source identifier is output in numerical form, where a value of 1 represents the component disturbance state, and a value of 0 represents the hardware drift state. This disturbance source identifier serves as the basis for determining which correction path to execute in subsequent steps.

[0052] Please see Figure 3 As shown, in S4, when the disturbance source indicator indicates a component disturbance, the current values ​​of the sound path length and timing reference remain unchanged, and only the deviation between the measured sound velocity and the theoretical sound velocity is converted into a compression factor fine-tuning amount according to a preset mapping relationship; when the indicator indicates hardware drift, the correction of the sound path length and timing reference is activated, while keeping the compression factor unchanged, and the corrected sound path length, timing reference, and current effective compression factor are output, specifically including: After obtaining the disturbance source identifier, the current disturbance type is determined based on the identifier's value. A value of 1 indicates a component disturbance state; a value of 0 indicates a hardware drift state. The corresponding correction processes are then performed according to each of these two states.

[0053] When the disturbance source identifier indicates a component disturbance state, calculate the weighted average of the measured sound velocities of each channel.

[0054] The weight value corresponding to each channel is preset according to the radial position of the channel in the pipe cross section. Specifically, the weight value of the channel located in the central area of ​​the pipe is greater than the weight value of the channel near the pipe wall, and the sum of all weight values ​​is equal to 1.

[0055] The measured sound velocity of each channel is multiplied by its corresponding weight value, and the results of all channel multiplications are summed to obtain the weighted average of the measured sound velocities. This weighted average is used as the representative sound velocity value of the cross section, and the difference is calculated with the theoretical sound velocity value based on the AGA 8 state equation. That is, the sound velocity deviation is obtained by subtracting the theoretical sound velocity value from the representative sound velocity value of the cross section.

[0056] A positive value for the sound speed deviation indicates that the measured sound speed is higher than the theoretical sound speed, while a negative value indicates that the measured sound speed is lower than the theoretical sound speed.

[0057] After obtaining the sound velocity deviation, its absolute value is used for subsequent interpolation calculations. Multiple continuous numerical intervals are pre-constructed, each interval having a starting boundary value, a ending boundary value, a starting fine-tuning value, and a ending fine-tuning value. The units of the starting and ending boundary values ​​are meters per second, and the starting and ending fine-tuning values ​​are dimensionless compression factor adjustments.

[0058] Specifically, the starting boundary value of the first numerical interval is 0, the ending boundary value is 0.1, the starting fine-tuning value is 0, and the ending fine-tuning value is 0.0005; The starting boundary value of the second numerical interval is 0.1, the ending boundary value is 0.5, the starting fine-tuning value is 0.0005, and the ending fine-tuning value is 0.002. The starting boundary value of the third numerical interval is 0.5, the ending boundary value is 1.0, the starting fine-tuning value is 0.002, and the ending fine-tuning value is 0.005. The starting boundary value of the fourth numerical interval is 1.0, the ending boundary value is 2.0, the starting fine-tuning value is 0.005, and the ending fine-tuning value is 0.01. The starting boundary value of the fifth numerical interval is 2.0, the ending boundary value is 5.0, the starting fine-tuning value is 0.01, and the ending fine-tuning value is 0.02.

[0059] The absolute value of the currently obtained sound speed deviation is compared with the starting and ending boundary values ​​of each numerical interval in turn to locate the target numerical interval into which the absolute value falls.

[0060] The specific positioning method is as follows: when the absolute value is greater than or equal to the starting boundary value of a certain interval and less than the ending boundary value of that interval, the interval is determined as the target value interval; when the absolute value is greater than or equal to the ending boundary value of the maximum interval, the maximum value interval is determined as the target value interval. After determining the target value interval, the relative position proportion is calculated using the following mathematical formula: Wherein, represents the absolute value of the sound speed deviation, represents the initial boundary value of the target numerical range, represents the final boundary value of the target numerical range, and represents the relative position percentage. The relative position percentage is a dimensionless value, ranging from 0 to 1.

[0061] After obtaining the relative position proportion, the compression factor adjustment amount is output using the following mathematical formula based on the starting and ending fine-tuning values ​​of the target value interval: Wherein, represents the initial fine-tuning value of the target numerical range, represents the final fine-tuning value of the target numerical range, represents the relative position percentage, and represents the output compression factor fine-tuning amount. The compression factor fine-tuning amount is superimposed on the theoretical compression factor; that is, the theoretical compression factor value is added to the compression factor fine-tuning amount to obtain the current effective compression factor.

[0062] Under component perturbation conditions, the current values ​​of the sound path length and timing reference remain unchanged, and no correction operations are performed.

[0063] When the disturbance source indicator indicates a hardware drift state, calculate the average deviation between the measured sound velocity and the theoretical sound velocity for each channel.

[0064] Specifically, the difference between the measured sound velocity and the theoretical sound velocity for each channel is calculated separately. The sum of the differences for all channels is then divided by the total number of channels to obtain the average deviation. The average deviation is multiplied by a first proportionality coefficient, which is 0.001, in seconds per meter, to obtain the path length correction. The mean deviation is multiplied by a second proportionality factor to obtain the timing reference correction, which is 0.0005 in seconds. The path length correction is then added to the current path length, i.e., the current path length value is added to the path length correction to obtain the corrected path length. The timing reference correction is then added to the current timing reference, i.e., the current timing reference value is added to the timing reference correction to obtain the corrected timing reference.

[0065] In the case of hardware drift, the theoretical compression factor remains unchanged, and is directly used as the current effective compression factor.

[0066] The current effective compression factor, the corrected sound path length, and the corrected timing reference obtained in the above steps are used as input parameters for subsequent steps.

[0067] If the current state is component perturbation, the output path length and timing reference are the uncorrected current values; if the current state is hardware drift, the output path length and timing reference are the corrected updated values. The current effective compression factor is output in both states.

[0068] In S5, based on the corrected path length, timing reference, and current effective compression factor, the average velocity of the cross section is reconstructed by combining the average velocity of each channel line, and the operating volumetric flow rate is converted to the standard volumetric flow rate under standard reference conditions as the output measurement result, specifically including: After obtaining the corrected path length, the corrected timing reference, and the current effective compression factor, the cross-sectional average velocity reconstruction and standard reference volumetric flow rate conversion are performed by combining the original measured values ​​of the average velocity of each channel and the original measurement time difference between the forward and reverse flows of each channel.

[0069] Based on the corrected path length, the path deviation of the line average velocity corresponding to each channel is corrected. Then, the timing deviation is inverted and corrected based on the corrected timing reference for the forward and reverse flow time difference corresponding to the corrected flow velocity of each channel to obtain the final flow velocity of each channel. Each channel is assigned a preset radial position weighting coefficient, and the sum of the products of the final flow velocity of each channel and its radial position weighting coefficient is used as the cross-sectional average flow velocity. Using the current effective compressibility factor as the compressibility factor value under standard reference conditions, the operating condition volumetric flow rate is converted into the standard reference volumetric flow rate according to the ratio between the standard reference pressure and the standard reference temperature, and the result is output as the measurement result.

[0070] During the sound path deviation correction process, the corrected sound path length and the original sound path length at factory calibration are obtained. The corrected sound path length is divided by the original sound path length to obtain the sound path correction ratio. The original measured value of the line average flow velocity corresponding to each channel is multiplied by the sound path correction ratio to obtain the corrected flow velocity for each channel. This corrected flow velocity eliminates the flow velocity measurement deviation introduced by the transducer sound path length changes with temperature or slight variations in installation position.

[0071] In the timing deviation inversion correction process, the corrected timing reference and the standard timing reference are first obtained.

[0072] The standard timing reference is the reference value corresponding to the timing frequency of the time-to-digital converter circuit in its factory-calibrated state, and this value is pre-stored in a non-volatile memory unit. The timing deviation ratio is obtained by dividing the corrected timing reference by the standard timing reference.

[0073] Take the reciprocal of the timing deviation ratio to obtain the inverse proportionality coefficient. Multiply the original downstream and upstream measurement time differences for each channel by this inverse proportionality coefficient to obtain the corrected downstream and upstream time differences.

[0074] This inverse proportional operation eliminates the time measurement deviation introduced by the drift of the timing crystal oscillator frequency.

[0075] After obtaining the corrected downstream time difference and the corrected upstream time difference, the difference between the corrected downstream time difference and the corrected upstream time difference is taken as the effective time difference.

[0076] Multiply the effective time difference by the corrected path length, and then divide by the current measured sound velocity of the channel to obtain the final time difference compensation for that channel.

[0077] This compensation amount reflects the time offset caused by the coupling effect of timing reference correction and sound path length correction. The final time difference compensation amount is calculated using the following formula: ; in, Indicates the first The final time difference compensation amount for the audio channel, in seconds; Indicates the first The time difference after channel correction, in seconds; Indicates the first The time difference of backflow after channel correction, in seconds; This indicates the corrected sound path length, in meters. Indicates the first The current measured sound velocity of the audio channel, in meters per second.

[0078] The final time difference compensation is added to the corrected downstream time difference and the corrected upstream time difference, respectively. That is, the corrected downstream time difference is added to the final time difference compensation to obtain the superimposed downstream time difference; the corrected upstream time difference is added to the final time difference compensation to obtain the superimposed upstream time difference.

[0079] The superimposed downstream time difference and the superimposed upstream time difference are used as the basis for calculating the final velocity of each channel. The following formula is used to calculate the final velocity of each channel: ; in, Indicates the first The final flow rate of the vocal tract, measured in meters per second; This indicates the corrected sound path length, in meters. Indicates the first The time difference of the current after channel superposition compensation, in seconds; Indicates the first The time difference of backflow after channel superposition compensation, in seconds.

[0080] This formula is based on the ultrasonic time-difference method for velocity measurement. It uses the relationship between the propagation time difference with and against the current and the sound path length to calculate the flow velocity. The product term in the denominator is used to eliminate the influence of sound speed changes on the flow velocity calculation.

[0081] After obtaining the final flow velocity of each channel, a preset radial position weighting coefficient is assigned to each channel. This weighting coefficient is predetermined based on the Gauss-Jacobi numerical integration method, reflecting the proportion of the annular region area of ​​the pipe cross-section represented by each channel to the total cross-sectional area, and the sum of all weighting coefficients equals 1. The final flow velocity of each channel is multiplied by its corresponding radial position weighting coefficient, and the products of all channels are summed to obtain the cross-sectional average flow velocity. The cross-sectional average flow velocity is multiplied by the pipe cross-sectional area to obtain the volumetric flow rate under operating conditions.

[0082] When converting to standard reference conditions, the current effective compression factor is used as the compression factor value under standard reference conditions.

[0083] Obtain the standard reference pressure and standard reference temperature values, where the standard reference pressure is 101.325 kPa and the standard reference temperature is 293.15 Kelvin.

[0084] The pressure ratio is obtained by dividing the standard reference pressure by the current pipeline static pressure; the temperature ratio is obtained by dividing the current medium temperature by the standard reference temperature.

[0085] Multiply the operating volumetric flow rate by the pressure ratio, then by the temperature ratio, and finally by the ratio of the current effective compressibility factor to the standard compressibility factor to obtain the standard reference volumetric flow rate at the standard reference pressure and standard reference temperature. The standard compressibility factor is set to 1.

[0086] The standard reference volumetric flow rate is output as the measurement result to complete the online correction and measurement process.

[0087] The working principle of this invention is as follows: The propagation time of ultrasonic waves in both directions is collected for each channel; the line-average velocity and measured line-average sound velocity of each channel are calculated; and the mole fraction of natural gas components and pipeline temperature and pressure parameters are obtained. The measured sound velocities of each channel are arranged radially along the pipeline and subjected to variational mode decomposition to separate low-frequency large-scale gradient components and high-frequency narrowband random components. Based on the low-frequency large-scale gradient components, the radial monotonicity index is calculated. This monotonicity index is matched with the upstream chromatographic component change timestamps within a time window, and an identifier distinguishing between component disturbances and hardware drift is output. When the identifier indicates component disturbance, the path length and timing reference remain unchanged; only the deviation between the measured and theoretical sound velocities is converted into a compression factor fine-tuning amount. When the identifier indicates hardware drift, the path length and timing reference are corrected while maintaining the compression factor unchanged. The corrected path length, timing reference, and current effective compression factor are output. Based on the corrected path length, timing reference, and current effective compression factor, and combined with the line-average velocity of each channel, the cross-sectional average velocity is reconstructed, and the operating volumetric flow rate is converted to the standard volumetric flow rate under standard reference conditions and output.

[0088] The foregoing has provided a detailed description of one embodiment of the present invention, but this description is merely a preferred embodiment and should not be construed as limiting the scope of the invention. All equivalent variations and modifications made within the scope of the claims of this invention should still fall within the patent coverage of this invention.

Claims

1. An online correction and metering method for ultrasonic natural gas flow meters, characterized in that, Includes the following steps: S1, collect the propagation time of ultrasonic waves in both directions of flow for each channel, calculate the line average flow velocity and the line average measured sound velocity for each channel, and simultaneously obtain the mole fraction of natural gas components and pipeline temperature and pressure parameters; S2, the measured sound velocities of each channel are arranged radially along the pipe, and variational mode decomposition is performed on them to separate the low-frequency large-scale gradient component that characterizes the component stratification effect and the high-frequency narrowband random component that characterizes the thermal expansion and contraction of the transducer and the drift of the timing reference. S3 calculates the radial monotonicity index based on the low-frequency large-scale gradient components, and matches the monotonicity index with the time stamp of the upstream chromatographic component change to output a disturbance source identifier to distinguish between state component disturbance and motion hardware drift. S4, when the disturbance source indicator indicates a component disturbance, keep the current values ​​of the sound path length and timing reference unchanged, and only convert the deviation between the measured sound velocity and the theoretical sound velocity into a compression factor fine-tuning amount according to a preset mapping relationship; When the indication is hardware drift, activate the correction of sound path length and timing reference, keep the compression factor unchanged, and output the corrected sound path length, timing reference and current effective compression factor. S5, based on the corrected sound path length, timing reference and current effective compression factor, reconstructs the cross-sectional average flow velocity by combining the average flow velocity of each channel line, and converts the working condition volumetric flow rate into the standard volumetric flow rate under standard reference conditions, as the output of the measurement result.

2. The online correction and metering method for ultrasonic natural gas flow meters according to claim 1, characterized in that, S2 specifically includes: The measured sound velocities of each channel are used as a spatial sequence that is discretely distributed along the radial direction of the pipe. The Hilbert transform is then applied to the spatial sequence to obtain the analytical signals of each modal component. The center frequency and bandwidth of each modal component are updated iteratively, and in each iteration, the upper limit of the bandwidth of the low-frequency large-scale gradient component is constrained to be less than the first preset value and its radial second-order difference norm is constrained to be less than the second preset value. At the same time, the lower limit of the bandwidth of the high-frequency narrowband random component is constrained to be greater than the third preset value. The iteration terminates when the sum of the center frequency changes of all modal components in two adjacent iterations is less than the fourth preset value, and the separated low-frequency large-scale gradient components and high-frequency narrowband random components are output.

3. The online correction and metering method for ultrasonic natural gas flow meters according to claim 2, characterized in that, The process of obtaining analytic signals for each modal component after performing a Hilbert transform on the spatial sequence specifically includes: Using the measured sound velocity of each channel as the initial real part of the analytical signal, the measured sound velocity distribution sequence along the radial direction is extended to generate a radially extended sequence. The Hilbert transform is applied to the radially extended sequence to generate an imaginary part sequence orthogonal to the real part. The real part sequence and the imaginary part sequence together constitute the analytical signal of each modal component. Extract the amplitude envelope corresponding to the original radial position in the analytical signal, perform differential operation between adjacent channels on the amplitude envelope, and output the amplitude change between each adjacent channel; The amplitude changes between each adjacent channel are compared with the preset radial continuity tolerance. When all amplitude changes are less than the radial continuity tolerance, the current analytical signal is output to the iterative update stage to participate in subsequent iterative calculations. Otherwise, radial smoothing is applied to the imaginary part sequence of the current analytical signal until the smoothed amplitude changes meet the radial continuity tolerance before outputting.

4. The online correction and metering method for ultrasonic natural gas flow meters according to claim 1, characterized in that, S3 specifically includes: For the low-frequency large-scale gradient component, the sound velocity difference between adjacent channels is calculated sequentially along the radial direction. The sign value of each sound velocity difference is accumulated, and the ratio of the absolute value of the accumulated result to the total number of channels is used as the monotonicity index output. Extract the moment when the mole fraction of a component jumps from the upstream chromatographic data as the component change timestamp, calculate the difference between the current time and the component change timestamp, and mark the time window as successfully matched when the absolute value of the difference is less than the preset time window value. When the monotonicity exponent is greater than the preset monotonicity threshold and the time window is successfully matched, the disturbance source identifier is set to the component disturbance state; otherwise, it is set to the hardware drift state, and the disturbance source identifier is output.

5. The online correction and metering method for ultrasonic natural gas flow meters according to claim 1, characterized in that, S4 specifically includes: The type of disturbance is determined based on the source identifier. When the identifier is a component disturbance, the measured sound velocities of each channel are weighted and averaged radially, and the difference is calculated with the theoretical sound velocity to obtain the sound velocity deviation. The corresponding compression factor fine-tuning amount is determined by piecewise linear interpolation based on the absolute value of the sound velocity deviation. The compression factor fine-tuning amount is then superimposed on the theoretical compression factor to output the current effective compression factor, while keeping the sound path length and timing reference at their current values. When hardware drift is identified, the average deviation between the measured sound velocity and the theoretical sound velocity of each channel is used as the correction basis. The average deviation is proportionally calculated and used as the path length correction amount and the timing reference correction amount, respectively. These are then superimposed on the current path length and the current timing reference, and the corrected path length and the corrected timing reference are output. At the same time, the theoretical compression factor is kept as the current effective compression factor. The current effective compression factor, the corrected path length, and the corrected timing reference are all used as input parameters for S5.

6. The online correction and metering method for ultrasonic natural gas flow meters according to claim 5, characterized in that, The determination of the corresponding compression factor fine-tuning amount after piecewise linear interpolation specifically includes: Multiple continuous numerical intervals between the absolute value of the sound velocity deviation and the compression factor fine-tuning amount are pre-constructed, and corresponding interval start fine-tuning value and interval end fine-tuning value are configured for each numerical interval. The absolute value of the current sound speed deviation is compared with the boundary value of each numerical interval to locate the target numerical interval into which the absolute value falls, and the relative position percentage of the absolute value within the target numerical interval is calculated. Based on the relative position ratio, interpolate linearly between the starting fine-tuning value and the ending fine-tuning value of the target value interval, output the compression factor fine-tuning amount corresponding to the absolute value, and superimpose the compression factor fine-tuning amount to the theoretical compression factor to obtain the current effective compression factor.

7. The online correction and metering method for ultrasonic natural gas flow meters according to claim 1, characterized in that, The conversion of the operating condition volumetric flow rate into the standard volumetric flow rate under standard reference conditions, and outputting it as the measurement result, specifically includes: Based on the corrected path length, the path deviation of the line average velocity corresponding to each channel is corrected to obtain the corrected velocity of each channel. Then, the timing deviation is inverted and corrected based on the time difference between the forward and reverse flow corresponding to the corrected velocity of each channel using the corrected timing reference to obtain the final velocity of each channel. Each channel is assigned a preset radial position weighting coefficient, and the sum of the products of the final flow velocity of each channel and its radial position weighting coefficient is used as the cross-sectional average flow velocity. Using the current effective compressibility factor as the compressibility factor value under standard reference conditions, the operating condition volumetric flow rate is converted to the standard reference volumetric flow rate at the standard reference pressure and standard reference temperature according to the pressure ratio and temperature ratio, and the standard reference volumetric flow rate is output as the measurement result.

8. The online correction and metering method for ultrasonic natural gas flow meters according to claim 7, characterized in that, The process of obtaining the final flow rate of each channel specifically includes: Obtain the deviation ratio between the corrected timing reference and the standard timing reference. Perform an inverse proportional operation on the original forward and reverse flow time differences corresponding to each channel with the deviation ratio to obtain the corrected forward flow time difference and the corrected reverse flow time difference. The difference between the corrected downstream time difference and the corrected upstream time difference is taken as the effective time difference. The effective time difference is multiplied by the corrected sound path length and then divided by the current measured sound speed to obtain the final time difference compensation amount of the sound channel. The final time difference compensation is superimposed on the corrected downstream time difference and the corrected upstream time difference. The superimposed downstream time difference and upstream time difference are used as the basis for calculating the final flow rate of each channel, and the final flow rate of each channel is output.

Images