RU2707576C1 - Method for calculating current difference of phase and frequency of signals of inertial flow meters (versions) - Google Patents

Method for calculating current difference of phase and frequency of signals of inertial flow meters (versions) Download PDF

Info

Publication number
RU2707576C1
RU2707576C1 RU2019113186A RU2019113186A RU2707576C1 RU 2707576 C1 RU2707576 C1 RU 2707576C1 RU 2019113186 A RU2019113186 A RU 2019113186A RU 2019113186 A RU2019113186 A RU 2019113186A RU 2707576 C1 RU2707576 C1 RU 2707576C1
Authority
RU
Russia
Prior art keywords
signal
parameters
frequency
signals
poles
Prior art date
Application number
RU2019113186A
Other languages
Russian (ru)
Inventor
Александр Сергеевич Семенов
Ольга Леонидовна Ибряева
Original Assignee
Федеральное государственное автономное образовательное учреждение высшего образования "Южно-Уральский государственный университет (национальный исследовательский университет)" ФГАОУ ВО "ЮУрГУ (НИУ)"
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Федеральное государственное автономное образовательное учреждение высшего образования "Южно-Уральский государственный университет (национальный исследовательский университет)" ФГАОУ ВО "ЮУрГУ (НИУ)" filed Critical Федеральное государственное автономное образовательное учреждение высшего образования "Южно-Уральский государственный университет (национальный исследовательский университет)" ФГАОУ ВО "ЮУрГУ (НИУ)"
Priority to RU2019113186A priority Critical patent/RU2707576C1/en
Application granted granted Critical
Publication of RU2707576C1 publication Critical patent/RU2707576C1/en

Links

Images

Classifications

    • GPHYSICS
    • G01MEASURING; TESTING
    • G01FMEASURING VOLUME, VOLUME FLOW, MASS FLOW OR LIQUID LEVEL; METERING BY VOLUME
    • G01F1/00Measuring the volume flow or mass flow of fluid or fluent solid material wherein the fluid passes through a meter in a continuous flow
    • G01F1/76Devices for measuring mass flow of a fluid or a fluent solid material
    • G01F1/78Direct mass flowmeters
    • G01F1/80Direct mass flowmeters operating by measuring pressure, force, momentum, or frequency of a fluid flow to which a rotational movement has been imparted
    • G01F1/84Coriolis or gyroscopic mass flowmeters

Abstract

FIELD: monitoring and measuring equipment.
SUBSTANCE: invention relates to control and measurement equipment and methods of processing one or more sensor signals in a flow meter and can be used in instrument-making during development and manufacturing of inertial flow meters. Method consists in using the signal representation from the inertial flow meter position sensors in the form of a sum of complex exponentials with unknown parameters. This allows after samples of signals to obtain a set of readings, each of which, in turn, will also be a sum of corresponding partial components with the same parameters, which can be presented in form of
Figure 00000056
, n=0…S>N, where lower digital index denotes sensor number, R1(2),k, ω1(2),k, α1(2),k – complex amplitude, frequency and attenuation of the k-th harmonic of the corresponding sensor,
Figure 00000057
M is the number of exponential components in the analyzed signal, N is the discrete number in the evaluation window, ε1(2),k(t) is the additive noise component of the corresponding signal. Used signal model enables to calculate the parameters R1(2),k, ω1(2),k, α1(2),k by method of matrix beams (method 1) or vector method of matrix beams (method 2), and from them to calculate required frequency and phase difference in accordance with expressions
Figure 00000058
For the second method, by the assumption of the poles identity, equations α1,k,n2,k,n and ω1,k,n2,k,n. Further, these parameters are transmitted to a flow parameter estimation system, after which the next n+1th unit of readings arrives into the system, which is equivalent to shifting the evaluation window by one sampling cycle, and the cycle of calculations is repeated again.
EFFECT: high accuracy of measuring mass and volume flow rate of a liquid medium in the presence of a perturbing phase (gas or solid) in inertial flow meters by reducing the number of samples and, accordingly, the time required to obtain the current flow rate.
2 cl, 9 dwg

Description

Изобретение относится к контрольно-измерительной технике и способам обработки одного или более сигналов датчиков в расходомере и может быть использовано в приборостроении при разработке и изготовлении кориолисовых расходомеров.The invention relates to instrumentation and methods for processing one or more sensor signals in a flow meter and can be used in instrumentation in the development and manufacture of Coriolis flow meters.

Массовые расходомеры Кориолиса широко используются для измерения массового расхода плотности и объемного расхода, а также получения другой информации о веществах, протекающим через трубопровод, как раскрыто в патенте US 4491025, МПК G01F 1/76, G01F 1/84, опубл. 01.01.1985 г. Эти расходомеры имеют одну или более расходомерных труб различных конфигураций. Каждая конфигурация трубы имеет набор мод собственных колебаний, в том числе, но не только, плоский изгиб, крутильную, радиальную и связанную моду. В типичном варианте применения измерения массового расхода методом Кориолиса используемая конфигурация труб возбуждается на одной или более колебательных мод при протекании вещества через трубопровод, и движение расходомерных труб измеряется в точках, разнесенных по длине трубы, как показано на фиг. 1Coriolis mass flowmeters are widely used to measure the mass flow rate of density and volumetric flow rate, as well as to obtain other information about substances flowing through the pipeline, as disclosed in patent US 4491025, IPC G01F 1/76, G01F 1/84, publ. 01/01/1985, These flowmeters have one or more flow tubes of various configurations. Each pipe configuration has a set of modes of natural vibrations, including, but not limited to, plane bending, torsional, radial, and coupled modes. In a typical Coriolis mass flow measurement application, the pipe configuration used is excited at one or more vibrational modes as the material flows through the pipeline, and the flow tubes are measured at points spaced along the length of the pipe, as shown in FIG. one

Приводной механизм возбуждает колебания расходомерной трубки. Когда нет вещества, протекающего через расходомер, все точки вдоль расходомерной трубки колеблются с идентичной фазой. По мере того как вещество начинает протекать через расходомерную трубку, ускорения Кориолиса приводят к тому, что каждая точка вдоль расходомерной трубки имеет различную фазу относительно других точек вдоль расходомерной трубки. Фаза на входной стороне расходомерной трубки запаздывает от приводного механизма, тогда как фаза на выходной стороне опережает приводной механизм.The drive mechanism excites vibrations of the flow tube. When there is no substance flowing through the flowmeter, all points along the flowmeter tube oscillate with an identical phase. As a substance begins to flow through the flow tube, Coriolis accelerations cause each point along the flow tube to have a different phase relative to other points along the flow tube. The phase on the input side of the flow tube is delayed from the drive mechanism, while the phase on the output side is ahead of the drive mechanism.

Датчики помещаются в различных точках на расходомерной трубке и преобразуют движение этих точек в соответствующий набор сигналов. Разность фаз соответствующих мод сигналов датчиков сигналов пропорциональна массовому расходу вещества, протекающего через расходомерную трубку или расходомерные трубки.The sensors are placed at various points on the flow tube and convert the movement of these points into an appropriate set of signals. The phase difference of the respective signal sensor signal modes is proportional to the mass flow rate of the substance flowing through the flow tube or flow tubes.

Из предыдущего уровня техники известны способы обработки сигналов для измерения массового расхода Кориолиса, основные из которых рассмотрены в [1]. Для оценки разности фаз необходимой моды колебаний используются дискретное преобразование Фурье (DFT), быстрое преобразование Фурье (FFT), различные модификации режекторных и следящих фильтров. Далее полученные оценки разности фаз и частоты колебаний блока расходомерных трубок могут быть использованы для того, чтобы вычислить массовый расход и плотность расходуемого вещества.Methods of signal processing for measuring Coriolis mass flow rate are known from the prior art, the main of which are considered in [1]. To estimate the phase difference of the required oscillation mode, discrete Fourier transform (DFT), fast Fourier transform (FFT), various modifications of notch and servo filters are used. Further, the obtained estimates of the phase difference and the oscillation frequency of the flow tube block can be used to calculate the mass flow rate and the density of the expendable substance.

Указанные способы работают достаточно хорошо в стационарном режиме, когда расходуемое вещество в расходомере является однородным и значения мод колебаний, определяемые плотностью и расходом контролируемого вещества, стабильны.These methods work quite well in a stationary mode, when the consumable in the flowmeter is homogeneous and the values of the vibration modes determined by the density and flow rate of the controlled substance are stable.

Однако, когда расходуемое вещество является неоднородным, например, в двухфазных потоках, где расходуемое вещество содержит жидкость и твердое вещество или имеются пузырьки газа в жидком расходуемом веществе, возникают быстрые флюктуации частот мод колебаний, которые не могут отслеживаться техникой предыдущего уровня, в силу принципиальных ограничений преобразования Фурье и основанных на нем фильтров.However, when the consumable is inhomogeneous, for example, in two-phase flows, where the consumable contains liquid and solid or there are gas bubbles in the liquid consumable, rapid fluctuations of the frequencies of the vibration modes occur that cannot be monitored by the technology of the previous level, due to fundamental limitations Fourier transforms and filters based on it.

Это эффект также является проблемой в режимах «пустой - полный - пустой», когда возможно многократное появление значительных объемов газа эквивалентных временному полному исчезновению измеряемой среды из расходомерных трубок.This effect is also a problem in the "empty - full - empty" modes, when the significant occurrence of significant volumes of gas is equivalent to the temporary complete disappearance of the measured medium from the flow tubes.

Описанные режимы приводят к некорректным измерениям частоты, что, в свою очередь влечет значительные погрешности (20% и более) в измерении разности фаз, а значит, в дальнейшем, плотности и расхода.The described modes lead to incorrect frequency measurements, which, in turn, entails significant errors (20% or more) in measuring the phase difference, and hence, in the future, density and flow rate.

Наиболее близким к заявляемому решению является способ, предложенный для решения указанной проблемы в патенте RU 2371678, МПК G01F 1/84 опубл. 27.10.2009 «Высокоскоростная оценка частоты и фазы расходомеров» и заключающийся в использовании преобразования Гильберта, реализованном на основе 90 - градусных фазовращателях. Способ обработки сигналов датчиков в расходомере для последующего вычисления массового расхода плотности и объемного расхода реализуется электронным измерительным оборудованием, содержащим интерфейс для приема первого сигнала датчика и второго сигнала датчика и связанную с интерфейсом систему обработки, предназначенную для формирования девяностоградусного сдвига фаз из первого сигнала датчика с помощью преобразования Гилберта и вычисления разности фаз из девяностоградусного сдвига фаз, первого сигнала датчика и второго сигнала датчика. Частоту вычисляют из первого сигнала датчика и девяностоградусного сдвига фаз. Второй девяностоградусный сдвиг фаз может быть сформирован из второго сигнала датчика.Closest to the claimed solution is the method proposed to solve this problem in the patent RU 2371678, IPC G01F 1/84 publ. 10.27.2009 "High-speed estimation of the frequency and phase of flowmeters" and consisting in the use of the Hilbert transform, implemented on the basis of 90-degree phase shifters. A method of processing sensor signals in a flow meter for subsequent calculation of the mass flow rate of density and volumetric flow rate is implemented by electronic measuring equipment comprising an interface for receiving a first sensor signal and a second sensor signal and an associated processing system for generating a ninety-degree phase shift from the first sensor signal using Hilbert transforms and calculates the phase difference from a ninety degree phase shift, the first sensor signal and the second signal Occupancy. The frequency is calculated from the first sensor signal and the ninety degree phase shift. A second ninety degree phase shift may be formed from a second sensor signal.

Способ позволяет существенно повысить динамические характеристики расходомера и уменьшить погрешность оценки разности фаз в условиях небольших и достаточно медленных вариации частоты.The method can significantly increase the dynamic characteristics of the flow meter and reduce the error in estimating the phase difference under conditions of small and rather slow frequency variations.

Тем не менее, при существенных и быстрых изменениях частоты колебаний расходомерных трубок, что характерно для режимов «полный - пустой - полный» или скачков плотности расходуемого вещества (резкое изменение газовой фазы в двухфазных средах), вышеперечисленные недостатки методов предшествующего уровня обнаруживает и указанный способ. Это связано с тем, что реализация широкополосных фазовращателей базируется на тех же методах обработки сигналов, что и методы предыдущего уровня техники, а значит и имеет те же принципиальные ограничения по динамике измерений, что и они. Кроме того, преобразование Гильберта ориентировано на работу с одной модой колебаний, что не выполняется при многофазном потоке. Для уменьшения влияния шумовых компонент в приведенном способе предлагается использовать фильтрацию, что дополнительно ухудшает быстродействие алгоритма и снижает эффективность его использования.However, with significant and rapid changes in the frequency of oscillations of the flow tubes, which is typical for the “full - empty - full” modes or jumps in the density of the consumed substance (a sharp change in the gas phase in two-phase media), the aforementioned disadvantages of the methods of the previous level are also detected by this method. This is due to the fact that the implementation of broadband phase shifters is based on the same signal processing methods as the methods of the prior art, and therefore has the same fundamental limitations on the dynamics of measurements as they are. In addition, the Hilbert transform is oriented to work with one mode of oscillation, which is not performed with multiphase flow. To reduce the influence of noise components in the above method, it is proposed to use filtering, which further affects the performance of the algorithm and reduces the efficiency of its use.

Несмотря на различные способы оценки разности фаз и частоты, приводящие к улучшению отдельных характеристик кориолисовых расходомеров, работающих с многофазными средами, все они, в том числе и использованный в качестве прототипа, имеют один и тот же недостаток. Эти преобразования можно отнести к непараметрическим методам [2], имеющими принципиальное ограничение на разрешение частот, связанное с временем наблюдения примерным соотношением Δω≈1/ΔT, где Δω и ΔT - необходимое разрешения по частоте и время наблюдения необходимое для его обеспечения, соответственно. Это соотношение накладывает жесткие требования на длительность наблюдаемого участка при требованиях повышенного разрешения, что в свою очередь ухудшает динамические характеристики алгоритмов обработки и затрудняет работу с нестационарными сигналами.Despite various methods for assessing the phase difference and frequency, leading to the improvement of individual characteristics of Coriolis flowmeters working with multiphase media, all of them, including those used as a prototype, have the same drawback. These transformations can be attributed to nonparametric methods [2], which have a fundamental restriction on the frequency resolution associated with the observation time by the approximate ratio Δω≈1 / ΔT, where Δω and ΔT are the necessary frequency resolutions and the observation time necessary to ensure it, respectively. This ratio imposes stringent requirements on the duration of the observed area with high resolution requirements, which in turn degrades the dynamic characteristics of processing algorithms and makes it difficult to work with non-stationary signals.

Задача, на решение которой в первую очередь направлено заявляемое изобретение, заключается в уменьшении погрешности измерения массового и объемного расхода жидкой среды при наличии возмущающей фазы (газовой или твердой) в кориолисовых расходомерах за счет сокращения количества отсчетов и, соответственно, времени, требуемого для получения текущего значения расхода.The problem to which the invention is primarily aimed is to reduce the error in measuring the mass and volumetric flow rate of a liquid medium in the presence of a disturbing phase (gas or solid) in Coriolis flow meters by reducing the number of readings and, accordingly, the time required to obtain the current flow rates.

Указанная задача достигается тем, что согласно первому варианту сигналы от датчиков положения расходомерных трубок кориолисова расходомера представляются в виде суммы комплексных экспонентThis task is achieved by the fact that according to the first embodiment, the signals from the position sensors of the flow tubes of the Coriolis flow meter are represented as the sum of complex exponents

Figure 00000001
Figure 00000001

где нижний цифровой индекс обозначает номер датчика, R1(2),k, ω1(2),k, α1(2),k - комплексная амплитуда, частота и затухание k-ой гармоники соответствующего датчика,

Figure 00000002
- полюса сигналов, М - число экспоненциальных компонент в исследуемом сигнале, ε1(2),k(t) - аддитивная шумовая компонента соответствующего сигнала. Для оценки неизвестных параметров до начала вычислений, на основе известного диапазона частот колебаний и динамических характеристик измерительной системы, определяют исходные параметры: период дискретизации сигналов Т, число отсчетов N, определяющее длительность текущего окна обработки и некоторый целочисленный параметр L,
Figure 00000003
и сохраняют их в системе обработки данных. Входные сигналы с 1-го и 2-го датчиков расходомера подают на указанную систему обработки, где сигналы дискретизируют с заданным периодом Т, и преобразуют в дискретную последовательность отсчетов от 1-го и 2-го датчиковwhere the lower digital index denotes the sensor number, R 1 (2), k , ω 1 (2), k , α 1 (2), k is the complex amplitude, frequency and attenuation of the kth harmonic of the corresponding sensor,
Figure 00000002
are the signal poles, M is the number of exponential components in the signal under study, ε 1 (2), k (t) is the additive noise component of the corresponding signal. To evaluate the unknown parameters before the start of calculations, based on the known range of oscillation frequencies and dynamic characteristics of the measuring system, the initial parameters are determined: the signal sampling period T, the number of samples N, which determines the duration of the current processing window and some integer parameter L,
Figure 00000003
and store them in a data processing system. Input signals from the 1st and 2nd sensors of the flow meter are fed to the specified processing system, where the signals are sampled with a given period T, and converted into a discrete sequence of samples from the 1st and 2nd sensors

Figure 00000004
Figure 00000004

n=0…S>N, n - номер текущего отсчета, которые используют для определения параметров сигнала первого и второго датчиков R1(2),k, z1(2),k, М, для каждого k=1, …, М и текущего отсчета n>N. Для вычисления параметров Rk, zk, М используют метод матричных пучков (ММП) [3], после чего вычисляют искомые параметры сигнала по формуламn = 0 ... S> N, n is the number of the current reference, which is used to determine the signal parameters of the first and second sensors R 1 (2), k , z 1 (2), k , M, for each k = 1, ..., M and the current reference n> N. To calculate the parameters R k , z k , M using the method of matrix beams (IMF) [3], and then calculate the desired signal parameters by the formulas

Figure 00000005
Figure 00000005

Далее определяют текущую разность фаз и частоту мод колебаний, необходимую для определения параметров потока, протекающего через расходомер, обычно массового расхода и плотности жидкости, соответствующие номеру n текущего отсчета. При этом частота моды и разность фаз по первому способу определяется по формулам (4)Next, determine the current phase difference and the frequency of the vibration modes, necessary to determine the parameters of the flow flowing through the flow meter, usually the mass flow rate and fluid density, corresponding to the number n of the current reference. Moreover, the frequency of the mode and the phase difference in the first method is determined by the formulas (4)

Figure 00000006
Figure 00000006

После этого, окно оценки параметра сдвигается на один дискрет, т.е. n1=n+1 и в качестве нового окна оценки используется новый массив отсчетов n1-N÷n1=n+1-N÷n+1 и вышеприведенная последовательность действий повторяется.After that, the parameter estimation window is shifted by one discrete, i.e. n1 = n + 1 and as a new evaluation window a new array of samples n1-N ÷ n1 = n + 1-N ÷ n + 1 is used and the above sequence of actions is repeated.

Согласно второму варианту, сигнал y1(t) и y2(t) от 1-го и 2-го датчиков расходомера, представляется суммой комплексных экспонентAccording to the second option, the signal y 1 (t) and y 2 (t) from the 1st and 2nd sensors of the flow meter is represented by the sum of complex exponentials

Figure 00000007
Figure 00000007

где нижний цифровой индекс обозначает номер датчика, R1(2),k - комплексная амплитуда, k-ой гармоники соответствующего датчика, а ωk, αk - частота и затухание, определяемые полюсом zk k-ой гармоники, одинаковым для обоих датчиков, М - число экспоненциальных компонент в исследуемом сигнале, ε1(2),k(t) - аддитивная шумовая компонента соответствующего сигнала, согласно изобретению изначально предполагается идентичность полюсов сигналов, т.е. z1k=z2k=zk, k=1..М, при этом, для расчета полюсов сигнала используются объединенные матрицы, получение в соответствии с векторным методом матричных пучков, что приводит к получению для каждого текущего отсчета n>N идентичных значений полюсов zk,n, k=1..М, которые используйся для расчета соответствующих комплексных амплитуд R1(2),k,n, после чего параметры сигнала определяются по формуламwhere the lower digital index denotes the sensor number, R 1 (2), k is the complex amplitude of the kth harmonic of the corresponding sensor, and ω k , α k is the frequency and attenuation determined by the pole z k of the kth harmonic, the same for both sensors , M is the number of exponential components in the signal under investigation, ε 1 (2), k (t) is the additive noise component of the corresponding signal, according to the invention, the identity of the signal poles is initially assumed, i.e. z 1k = z 2k = z k , k = 1..M, in this case, the combined matrices are used to calculate the signal poles, obtaining matrix beams in accordance with the vector method, which leads to obtaining identical pole values for each current reference n> N z k, n , k = 1.. M, which are used to calculate the corresponding complex amplitudes R 1 (2), k, n , after which the signal parameters are determined by the formulas

Figure 00000008
Figure 00000008

и далее значения ωk,n и ϕk,n1,k,n2,k,n записываются в систему обработки данных оценки параметров потока, после чего в систему поступает очередной, n+1 - й блок отсчетов, что эквивалентно сдвигу окна оценки на один такт дискретизации, и цикл вычислений повторяется снова.and further, the values of ω k, n and ϕ k, n = ϕ 1, k, n2, k, n are written to the data processing system for estimating the flow parameters, after which the next, n + 1-th block of samples arrives at the system, which is equivalent to shifting the estimation window by one sampling clock, and the calculation cycle is repeated again.

Уменьшение времени обработки происходит за счет того, что в параметрических методах минимальное число отсчетов определяется не требуемым разрешением по частоте, как в непараметрических методах, а числом неизвестных параметров, которое может быть в десятки и сотни раз меньше, чем число отсчетов на интервале разрешения частоты методами Фурье и другими подобными методам. При этом, в первом приближении, длительность дискреты в параметрическом методе не регламентируется, т.е. формально необходимое число отсчетов может быть получено за интервал меньший, чем период колебании самой высокочастотной моды, в то время как для непараметрических методов необходимый интервал наблюдения может составлять десятки периодов. Учитывая это, предлагаемые способы выдают некоторую «точечную» оценку параметров в окне (на интервале наблюдения), получение зависимости этих сигналов от времени реализуется путем смещения окна (интервала наблюдения) при получении нового отсчета. Алгоритм смещения окна при получении нового отсчета отображают соответствующие шаги, показанные на фиг. 2 и 3.The reduction of the processing time is due to the fact that in parametric methods the minimum number of samples is determined not by the required frequency resolution, as in non-parametric methods, but by the number of unknown parameters, which can be tens or hundreds of times less than the number of samples on the frequency resolution interval by methods Fourier and other similar methods. At the same time, as a first approximation, the duration of the discretes in the parametric method is not regulated, i.e. Formally, the required number of samples can be obtained for an interval shorter than the period of oscillation of the highest frequency mode, while for nonparametric methods, the required observation interval can be tens of periods. Given this, the proposed methods give some "point" estimate of the parameters in the window (on the observation interval), obtaining the dependence of these signals on time is realized by shifting the window (observation interval) when receiving a new reference. The window shift algorithm upon receipt of a new sample is represented by the corresponding steps shown in FIG. 2 and 3.

Сущность изобретения поясняется следующими графическими материалами:The invention is illustrated by the following graphic materials:

Фиг. 1 Блок-схема кориолисова расходомера, где позициями обозначены следующие элементы: корпус расходомера 1, расходомерные трубки 2, датчики 3;FIG. 1 Block diagram of a Coriolis flowmeter, where the following elements are indicated by positions: flowmeter housing 1, flow tubes 2, sensors 3;

Фиг. 2 Блок-схема алгоритма по 1 варианту способа;FIG. 2 Block diagram of the algorithm according to 1 variant of the method;

Фиг. 3 Блок-схема алгоритма по 2 варианту способа;FIG. 3 Block diagram of the algorithm according to the 2nd variant of the method;

Фиг. 4 Сравнение алгоритмов в одномодовом режиме работы; отображены результаты расчетов по первому варианту способа и с использованием преобразования Гильберта по патенту RU 2371678, а также фактическое значение измеряемых параметров. Для демонстрации зависимости от длительности окна оценки графики и расчеты приведены для N=25 и N=100. В соответствующей таблице приведены оценки погрешности методов при измерении частоты, разности фаз и амплитуды сигналов.FIG. 4 Comparison of algorithms in single-mode operation; the results of calculations according to the first embodiment of the method and using the Hilbert transform according to the patent RU 2371678 are displayed, as well as the actual value of the measured parameters. To demonstrate the dependence on the window duration, the graph estimates and calculations are given for N = 25 and N = 100. The corresponding table provides estimates of the error of the methods in measuring the frequency, phase difference, and signal amplitude.

Фиг. 5 Сравнение алгоритмов при наличии второй моды; отображены результаты расчетов по второму варианту способа и с использованием преобразования Гильберта по патенту RU 2371678, а также фактическое значение измеряемых параметров при наличии дополнительной синусоидальной компоненты, являющейся в данном случае помехой, т.е моделируется двухмодовый режим работы. Добавленная синусоида имеет частоту 100 Гц. И амплитуду 20 dB от основной. Расчеты погрешностей методов для оцениваемых параметров приведены в таблице.FIG. 5 Comparison of algorithms in the presence of a second mode; the calculation results are shown according to the second variant of the method and using the Hilbert transform of the patent RU 2371678, as well as the actual value of the measured parameters in the presence of an additional sinusoidal component, which in this case is an obstacle, i.e., a two-mode operation mode is simulated. The added sine wave has a frequency of 100 Hz. And the amplitude is 20 dB from the main one. Calculation of error methods for the estimated parameters are given in the table.

Фиг. 6 Оценки и истинное значения полюсов мод, рассчитанные при помощи варианта 1 (верхняя и средняя диаграммы) и варианта 2 (нижняя диаграмма). Истинное значение полюсов z_1=0,8070+0,5863j и z_2=0,8484+0,5199j отмечено кружочком, а оцененное крестиком. Концентрация оценок вокруг истинных значений на нижних рисунках свидетельствуете об уменьшении погрешности при использовании второго варианта способа, что подтверждают данные приведенные в соответствующей таблице.FIG. 6 Estimates and true values of the mode poles calculated using option 1 (upper and middle diagrams) and option 2 (lower diagram). The true value of the poles z_1 = 0.8070 + 0.5863j and z_2 = 0.8484 + 0.5199j is marked with a circle, and estimated with a cross. The concentration of estimates around the true values in the lower figures indicates a decrease in the error when using the second variant of the method, which is confirmed by the data given in the corresponding table.

Фиг. 7 Частота сигналов, рассчитанная по варианту 1 (сверху и посередине) и варианту 2 (снизу). Отметим значительное уменьшение дисперсии оценки частоты на нижнем графике (вариант 2);FIG. 7 Signal frequency calculated according to option 1 (top and middle) and option 2 (bottom). We note a significant decrease in the variance of the frequency estimate in the lower graph (option 2);

Фиг. 8 Амплитуда, рассчитанная по варианту 1 (светлым) и варианту 2 (темным), для правой и левой катушки (сигнал 1, сигнал 2). Видно значительное уменьшение дисперсии оценки амплитуды при использовании второго варианта способа;FIG. 8 Amplitude calculated according to option 1 (light) and option 2 (dark) for the right and left coils (signal 1, signal 2). You can see a significant decrease in the variance of the amplitude estimate when using the second variant of the method;

Фиг. 9 Разница фаз сигналов, вычисленная по способу 1 (синим) и способу 2 (красным). Отметим уменьшение дисперсии оценки разности фаз, что свидетельствует об уменьшении погрешности при использовании способа 2.FIG. 9 The phase difference of the signals calculated by method 1 (blue) and method 2 (red). Note the decrease in the variance of the phase difference estimate, which indicates a decrease in the error when using method 2.

В качестве основы модели сигнала расходомера, работающего в условиях двухфазного потока использовалась модель, предложенная в [2].The model proposed in [2] was used as the basis of the signal model of a flowmeter operating in a two-phase flow.

При двухмодовом режиме работы мощность «мешающей» моды принималось 20 dB от основной.In the two-mode operation mode, the power of the “interfering” mode was assumed to be 20 dB from the main one.

Алгоритм заявляемого способа по первому варианту приведен на фиг. 2 и содержит следующие шаги:The algorithm of the proposed method according to the first embodiment is shown in FIG. 2 and contains the following steps:

1. Установка исходных параметров алгоритма, их сохранение в вычислительной системе, дискретизация сигналов и передача каждого отсчета по мере его получения в вычислительную систему.1. Setting the initial parameters of the algorithm, their storage in the computer system, discretization of signals and the transmission of each sample as it is received in the computer system.

2. Формирование массива («окна») из N>=4M последовательных отсчетов n-М÷n по каждому датчику, где n - текущий отсчет.2. Array formation (“windows”) from N> = 4M consecutive samples n-М ÷ n for each sensor, where n is the current sample.

3. Формирование матриц Y a , Yb (5), (6) для каждого датчика.3. Formation of matrices Y a , Y b (5), (6) for each sensor.

4. Нахождение SVD-разложения матриц Y a .4. Finding the SVD decomposition of the matrices Y a .

5. Оценивание число М полюсов сигнала.5. Estimation of the number M of signal poles.

6. Нахождение усеченных до ранга М псевдообратных матриц

Figure 00000009
(10).6. Finding truncated to rank M pseudoinverse matrices
Figure 00000009
(10).

7. Оценивание полюсов zk для каждого датчика путем вычисления собственных значений матриц ZE (12).7. Estimation of the poles z k for each sensor by calculating the eigenvalues of the matrices Z E (12).

8. Оценивание амплитуд Rk для каждого датчика, путем решения матричного уравнения (13) методом наименьших квадратов.8. Estimation of the amplitudes R k for each sensor by solving the matrix equation (13) by the least squares method.

9. Вычисление параметров мод в соответствии с (14), (15).9. Calculation of mode parameters in accordance with (14), (15).

10. Смещение «окна» на одну позицию, т.е. n=n+1.10. The shift of the "window" by one position, ie n = n + 1.

На первом шаге способа устанавливаются необходимый период дискретизации сигналов Т и размер «окна» обработки N>4M. После их сохранения в вычислительной системе с помощью АЦП реализуется процесс дискретизации и получения отсчетов сигнала, которые по мере их поступления также передаются в вычислительную систему.At the first step of the method, the necessary sampling period of the signals T and the size of the processing window N> 4M are set. After they are stored in the computer system using the ADC, the process of sampling and obtaining samples of the signal is implemented, which, as they arrive, are also transmitted to the computer system.

Далее, в соответствии с ММП, находят полюсы zk с помощью решения задачи на обобщенные собственные значения пучка матриц, составленных из отсчетов сигнала у(n). По найденным полюсам ММП находят оценки для Rk и параметры мод колебаний. Изложим подробнее способ с использованием ММП.Further, in accordance with the IMF, the poles z k are found by solving the problem of generalized eigenvalues of a matrix beam composed of signal samples y (n). Based on the found poles of the IMF, estimates for R k and vibration mode parameters are found. We describe in more detail the method using MMP.

После дискретизации сигналов с датчиков и формирования массивов измерений по каждому датчику (шаг 2), формируют две матрицы Y a , Yb (шаг 3) размеров (N-L)×L, следующим образом:After sampling the signals from the sensors and forming measurement arrays for each sensor (step 2), two matrices Y a , Y b (step 3) of dimensions (NL) × L are formed, as follows:

Figure 00000010
Figure 00000010

Здесь М≤L≤N-М параметр метода. Показано [3], что при выборе

Figure 00000011
дисперсия оценки полюсов zk будет минимальна, т.е. ММП будет наименее чувствителен к шуму. Для конкретизации можно положить
Figure 00000012
(функция Round округляет число до ближайшего целого).Here M≤L≤N-M parameter of the method. It is shown [3] that when choosing
Figure 00000011
the variance of the estimate of the poles z k will be minimal, i.e. MMP will be the least sensitive to noise. For concretization, we can put
Figure 00000012
(the Round function rounds a number to the nearest integer).

Основой метода матричных пучков является сингулярное разложение (SVD разложение).The basis of the matrix beam method is singular decomposition (SVD decomposition).

Для матриц Y a , Yb справедлива следующая факторизация [3]:For the matrices Y a , Y b the following factorization is true [3]:

Figure 00000013
Figure 00000013

ЗдесьHere

Figure 00000014
Figure 00000014

Figure 00000015
Figure 00000015

Рассмотрим пучок матриц Yb-λY a =ZLR(Z-λE)ZR, где Е - единичная матрица -го порядка. Можно показать [3], что при М≤L≤N-М и λ≠zk, k=1, 2, …, М, ранг матрицы Yb-λY a равен М. Однако, если λ=zk, то k-ая строка матрицы Z-λЕ - нулевая и ранг этой матрицы равен М-1.Consider a matrix pencil Y b -λY a = Z L R (Z-λE) Z R , where E is the ith order identity matrix. It can be shown [3] that for M≤L≤N-M and λ ≠ z k , k = 1, 2, ..., M, the rank of the matrix Y b -λY a is M. However, if λ = z k , then The k-th row of the matrix Z-λЕ is zero and the rank of this matrix is M-1.

Таким образом, полюсы

Figure 00000016
могут быть найдены как обобщенные собственные значения пучка матриц Yb-λY a [3], т.е. как собственные значения матрицы
Figure 00000017
(Здесь индекс
Figure 00000018
используется для обозначения псевдообратной матрицы.)So the poles
Figure 00000016
can be found as generalized eigenvalues of the matrix pencil Y b -λY a [3], ie as eigenvalues of a matrix
Figure 00000017
(Here is the index
Figure 00000018
used to denote a pseudoinverse matrix.)

В случае зашумленных данных матрицу Y a для эффективной фильтрации шума необходимо подвергнуть операции сингулярного (SVD) разложенияIn the case of noisy data, the matrix Y a must be subjected to singular (SVD) decomposition for effective noise filtering

Figure 00000019
Figure 00000019

Здесь U, V - унитарные матрицы, S - диагональная матрица, элементами которой являются сингулярные числа матрицы Y a . (Верхний индекс Т означает операцию транспонирования.)Here U, V are unitary matrices, S is the diagonal matrix, whose elements are singular numbers of the matrix Y a . (The superscript T stands for transpose.)

Заметим, что в случае отсутствия шума диагональная матрица S имеет ровно М ненулевых сингулярных чисел, все последующие равны нулю. В случае зашумленного сигнала ненулевых сингулярных чисел уже не будет, однако между первыми М и последующими сингулярными числами матрицы S будет наблюдаться ярко выраженный скачок, который и позволит определить число комплексных экспонент в сигнале.Note that in the absence of noise, the diagonal matrix S has exactly M nonzero singular numbers, all subsequent ones are equal to zero. In the case of a noisy signal, there will no longer be non-zero singular numbers, however, between the first M and subsequent singular numbers of the matrix S, a pronounced jump will be observed, which will determine the number of complex exponentials in the signal.

Сингулярное разложение матрицы Y a позволяет определить число истинных экспонент сигнала М. Кроме того, оно может быть использовано для нахождения псевдообратной матрицы

Figure 00000020
На практике используется усеченная до ранга М псевдообратная матрица:The singular decomposition of the matrix Y a allows us to determine the number of true exponents of the signal M. In addition, it can be used to find the pseudoinverse matrix
Figure 00000020
In practice, a pseudoinverse matrix truncated to rank M is used:

Figure 00000021
Figure 00000021

где σ1, …, σM - это М наибольших сингулярных чисел матрицы

Figure 00000022
соответствующие им сингулярные векторы,
Figure 00000023
, S0=diag(σ1, …, σM).where σ 1 , ..., σ M is the M largest singular numbers of the matrix
Figure 00000022
their corresponding singular vectors,
Figure 00000023
, S 0 = diag (σ 1 , ..., σ M ).

После нахождения матрицы

Figure 00000024
для оценки полюсов zk, k=1, …, М, сигнала остается только найти М собственных чисел матрицы
Figure 00000025
или, в силу следующей цепочки равенств:After finding the matrix
Figure 00000024
to estimate the poles z k , k = 1, ..., M, the signal can only be found M eigenvalues of the matrix
Figure 00000025
or, by virtue of the following chain of equalities:

Figure 00000026
Figure 00000026

найти М собственных чисел матрицыfind M eigenvalues of the matrix

Figure 00000027
Figure 00000027

Далее, в методе матричных пучков по известным М и zk находятся комплексные амплитуды Rk из решения следующей задачи:Further, in the method of matrix beams from the known M and z k are complex amplitudes R k from solving the following problem:

Figure 00000028
Figure 00000028

После вычисления комплексных амплитуд Rk можно определить искомые параметры сигнала по формуламAfter calculating the complex amplitudes R k, it is possible to determine the desired signal parameters by the formulas

Figure 00000029
Figure 00000029

Далее, на шаге 9 вычисляют текущую разность фаз и частоту мод колебаний, необходимую для определения параметров потока, протекающего через расходомер, обычно массового расхода и плотности жидкости, по формуламNext, in step 9, the current phase difference and the frequency of the oscillation modes, which are necessary to determine the parameters of the stream flowing through the flowmeter, usually the mass flow rate and the density of the liquid, are calculated by the formulas

Figure 00000030
Figure 00000030

На этом шаге вычисление параметров мод колебаний для отсчета n заканчивается, и они передаются в систему оценки параметров потока, после чего на шаге 10 в систему поступает очередной, n+1 - й блок отсчетов, что эквивалентно сдвигу окна оценки на один такт дискретизации.At this step, the calculation of the parameters of the vibrational modes for the sample n ends, and they are transferred to the system for estimating the flow parameters, after which at step 10 the next, n + 1-st block of samples enters the system, which is equivalent to shifting the estimation window by one sampling cycle.

Модификации метода и иные способы реализации ММП не изменяют сущности варианта.Modifications of the method and other methods for implementing IMF do not change the essence of the variant.

Согласно второму варианту способ реализуется путем использования информации об идентичности полюсов сигнала с датчиков, предложенным в [5].According to the second variant, the method is implemented by using information about the identity of the poles of the signal from the sensors proposed in [5].

Алгоритм второго варианта способа приведен на фиг. 3 и содержит следующие шаги:The algorithm of the second variant of the method is shown in FIG. 3 and contains the following steps:

1. Установка исходных параметров алгоритма, их сохранение в вычислительной системе, дискретизация сигналов и передача каждого отсчета по мере его получения в вычислительную систему.1. Setting the initial parameters of the algorithm, their storage in the computer system, discretization of signals and the transmission of each sample as it is received in the computer system.

2. Формирование массива («окна») из N>=4M последовательных отсчетов n-М÷n по каждому датчику, где n-текущий отсчет.2. Array formation (“windows”) from N> = 4M consecutive samples n-М ÷ n for each sensor, where n is the current sample.

3. Формирование матриц Y a 1, Yb1 Y a 2, Yb2 (5), (6) для каждого датчика.3. Formation of matrices Y a 1 , Y b1 Y a 2 , Y b2 (5), (6) for each sensor.

4. Формирование объединенных матриц пучка

Figure 00000031
4. Formation of combined beam matrices
Figure 00000031

5. Нахождение SVD-разложения матрицы Y a U.5. Finding the SVD decomposition of the matrix Y a U.

6. Оценивание числа М полюсов сигнала.6. Estimation of the number M of signal poles.

7. Нахождение усеченной до ранга М объединенной псевдообратной матрицы

Figure 00000032
(10).7. Finding a truncated to rank M joint pseudoinverse matrix
Figure 00000032
(10).

8. Оценивание М объединенных полюсов zk путем вычисления собственных значений матрицы ZE, полученной на основе (12) с использованием усеченной матрицы

Figure 00000033
8. Estimating M of the combined poles z k by calculating the eigenvalues of the matrix Z E obtained on the basis of (12) using the truncated matrix
Figure 00000033

9. Оценивание амплитуд Rk для каждого датчика, путем решения матричного уравнения (13) методом наименьших квадратов с подстановкой значений полюсов zk, полученных на предыдущем шаге.9. Estimation of the amplitudes R k for each sensor by solving the matrix equation (13) by the least square method with the substitution of the values of the poles z k obtained in the previous step.

10. Вычисление параметров мод в соответствии с (14), (15) (с учетом получения параметров частоты и затухания по объединенной матрице).10. Calculation of the mode parameters in accordance with (14), (15) (taking into account obtaining the frequency and attenuation parameters from the combined matrix).

11. Смещение «окна» на одну позицию, т.е. n=n+1.11. The shift of the "window" by one position, ie n = n + 1.

Основным отличием второго варианта от первого варианта способа является вычисление общих полюсов сигнала по объединенным матрицам пучка. В связи с этим первые три шага варианта совпадают с первыми тремя шагами основного способа.The main difference between the second variant and the first variant of the method is the calculation of the common signal poles from the combined beam matrices. In this regard, the first three steps of the variant coincide with the first three steps of the main method.

Далее, в целях получения оценок объединенных полюсов составляются блочные матрицы пучка

Figure 00000034
Как показано в [5], М полюсов zk являются обобщенными собственными значениями пучка матриц YbU-λY a U. В связи с этим значения объединенных полюсов могут быть получены в соответствии с алгоритмами первого варианта способа путем замены матриц Y a , Yb на матрицы Y a U, YbU. То есть шаги 5, 6, 7, 8 второго варианта аналогичны шагам 4, 5, 6, 7 первого варианта способа при соответствующей замене матриц. Комплексные амплитуды Rmk на шаге 9 алгоритма находятся затем из решения 2 систем уравнений:Further, in order to obtain estimates of the combined poles, block beam matrices are compiled
Figure 00000034
As shown in [5], the M poles z k are generalized eigenvalues of the matrix pencil Y bU −λY a U. In this regard, the values of the combined poles can be obtained in accordance with the algorithms of the first embodiment of the method by replacing the matrices Y a , Y b with the matrices Y a U , Y bU . That is, steps 5, 6, 7, 8 of the second embodiment are similar to steps 4, 5, 6, 7 of the first embodiment of the method with the corresponding replacement of the matrices. The complex amplitudes R mk at step 9 of the algorithm are then found from the solution of 2 systems of equations:

Figure 00000035
Figure 00000035

где z1…zM - объединенные полюса, полученные на шаге 8 алгоритма. Вычисление интересующих параметров мод по полученным значениям Rik и zk, i=1..M, k=1, 2 проводится на шаге 10 по выражениям (14), (15) аналогично первому варианту способа, с учетом идентичности полюсов для первого и второго сигнала.where z 1 ... z M are the combined poles obtained in step 8 of the algorithm. The calculation of the mode parameters of interest based on the obtained values of R ik and z k , i = 1..M, k = 1, 2, is carried out in step 10 using expressions (14), (15) similarly to the first method variant, taking into account the identity of the poles for the first and second signal.

Для сравнения характеристик предлагаемого алгоритма и прототипа было проведено математическое моделирование их работы. Для генерации исходных сигналов была использована модель, предложенная в [4] со следующими параметрами: основная частота изменяется в пределах 85-100 Гц, амплитуда основной частоты изменяется в пределах 0.05-0.3 В, разность фаз - в пределах 0-4 град., частота дискретизации 2 кГц. Для устранения резких изменений параметров их реализации были пропущены через низкочастный фильтр с полосой 6 Гц.To compare the characteristics of the proposed algorithm and the prototype, mathematical modeling of their work was carried out. To generate the initial signals, we used the model proposed in [4] with the following parameters: the fundamental frequency varies between 85-100 Hz, the amplitude of the fundamental frequency varies between 0.05-0.3 V, the phase difference is within 0-4 deg., Frequency 2 kHz sampling. To eliminate sudden changes in the parameters of their implementation, they were passed through a low-pass filter with a band of 6 Hz.

Выводы. Из приведенных результатов расчетов следует, что предложенный способ решает задачу по уменьшению погрешности измерения частоты и разности фаз при малых временах измерения, причем как показывает моделирование при уменьшении времени измерения преимущество предлагаемого метода увеличивается. Наличие дополнительных мод сокращает разницу между методами, но и в этом случае преимущество предлагаемого метода остается существенным. Уменьшение погрешности измерения частоты и разности фаз позволяют предположить, что соответствующим образом будет уменьшена и погрешность оценки расхода и плотности контролируемой жидкости при двухфазных режимах работы. Учет информации об идентичности полюсов согласно варианта 2 может существенно уменьшить погрешность измерения параметров. Однако при этом необходимо учитывать возможное усложнение системы обработки за счет увеличения размерности матриц, а также конкретные варианты конструкции расходомера, возможно приводящие к не идентичности потоков в расходомерных трубках.Conclusions. From the above calculation results, it follows that the proposed method solves the problem of reducing the error in measuring the frequency and phase difference at short measurement times, and as the simulation shows with decreasing measurement time, the advantage of the proposed method increases. The presence of additional modes reduces the difference between the methods, but in this case, the advantage of the proposed method remains significant. The decrease in the error in measuring the frequency and phase difference allows us to assume that the error in estimating the flow rate and density of the controlled liquid under two-phase operation modes will be correspondingly reduced. Accounting for information on the identity of the poles according to option 2 can significantly reduce the measurement error of the parameters. However, it is necessary to take into account the possible complication of the processing system by increasing the dimension of the matrices, as well as specific design options for the flowmeter, possibly leading to non-identical flows in the flow tubes.

Список литературыList of references

1. М. Li and М. Henry, "Signal processing methods for Coriolis Mass Flow Metering in two-phase flow conditions," in 2016 IEEE International Conference on Industrial Technology (ICIT), 2016, pp. 690-695.1. M. Li and M. Henry, "Signal processing methods for Coriolis Mass Flow Metering in two-phase flow conditions," in 2016 IEEE International Conference on Industrial Technology (ICIT), 2016, pp. 690-695.

2. Marple S.L. Digital spectral analysis: with applications. Prentice-Hall, 1987. 492 p.2. Marple S.L. Digital spectral analysis: with applications. Prentice-Hall, 1987.492 p.

3. Hua Y., Sarkar Т.K. Matrix Pencil Method for Estimating Parameters of Exponentially Damped/Undamped Sinusoids in Noise // IEEE Trans. Acoust. 1990. Vol. 38, №5. P. 814-824.3. Hua Y., Sarkar T.K. Matrix Pencil Method for Estimating Parameters of Exponentially Damped / Undamped Sinusoids in Noise // IEEE Trans. Acoust 1990. Vol. 38, No. 5. P. 814-824.

4. M. Li and M. Henry, "Complex Bandpass Filtering for Coriolis Mass Flow Meter Signal Processing," in Industrial Electronics Society (IECON), 2016, pp. 4952-4957.4. M. Li and M. Henry, "Complex Bandpass Filtering for Coriolis Mass Flow Meter Signal Processing," in Industrial Electronics Society (IECON), 2016, pp. 4952-4957.

5. M.P. Henry, O.L. Ibryaeva, D.D. Salov A.S.S. Matrix Pencil Method for Estimation of Parameters of Vector Processes // Bull. South Ural State Univ. Ser. "Mathematical Model. Program. Comput. Softw. 2017. Vol. 10, №4. P. 92-104.5. M.P. Henry, O.L. Ibryaeva, D.D. Salov A.S.S. Matrix Pencil Method for Estimation of Parameters of Vector Processes // Bull. South Ural State Univ. Ser. "Mathematical Model. Program. Comput. Softw. 2017. Vol. 10, No. 4. P. 92-104.

Claims (12)

1. Способ оценки текущей разности фаз и частоты сигналов y1(t) и y2(1) от 1-го и 2-го датчиков расходомера, представляемых суммой комплексных экспонент1. A method for estimating the current phase difference and frequency of signals y 1 (t) and y 2 (1) from the 1st and 2nd flowmeter sensors, represented by the sum of complex exponents
Figure 00000036
Figure 00000036
где нижний цифровой индекс обозначает номер датчика, R1(2),k, ω1(2),k, α1(2),k - комплексная амплитуда, частота и затухание k-й гармоники соответствующего датчика, полюса сигналов
Figure 00000037
М - число экспоненциальных компонент в исследуемом сигнале, ε1(2),k(t) - аддитивная шумовая компонента соответствующего сигнала, заключающийся в том, что предварительно на основе известного диапазона частот колебаний и динамических характеристик измеряемой системы определяют исходные параметры: период дискретизации сигналов Т, число отсчетов N, определяющее длительность текущего окна обработки и некоторый целочисленный параметр L,
Figure 00000038
которые сохраняют в системе обработки данных, после чего входные сигналы с 1-го и 2-го датчиков расходомера подают на систему обработки, где сигналы дискретизируют с заданным периодом Т и преобразуют в дискретную последовательность отсчетов от 1-го и 2-го датчиков
Figure 00000039
where the lower digital index denotes the sensor number, R 1 (2), k , ω 1 (2), k , α 1 (2), k is the complex amplitude, frequency and attenuation of the kth harmonic of the corresponding sensor, signal pole
Figure 00000037
M is the number of exponential components in the signal under investigation, ε 1 (2), k (t) is the additive noise component of the corresponding signal, which consists in the fact that preliminary parameters are determined based on the known range of vibration frequencies and dynamic characteristics of the measured system: the signal sampling period T, the number of samples N, which determines the duration of the current processing window and some integer parameter L,
Figure 00000038
which are stored in the data processing system, after which the input signals from the 1st and 2nd sensors of the flow meter are fed to the processing system, where the signals are sampled with a given period T and converted into a discrete sequence of samples from the 1st and 2nd sensors
Figure 00000039
n=0…S>N, n - номер текущего отсчета, массив которых используют для определения комплексных амплитуды, значения полюсов и числа полюсов R1(2),k, z1,(2),k, М, для каждого k=1, …, М и текущего отсчета n>N методом матричных пучков (ММП), после чего находят искомые параметры сигнала по формулам
Figure 00000040
n = 0 ... S> N, n is the number of the current reference, an array of which is used to determine the complex amplitudes, the values of the poles and the number of poles R 1 (2), k , z 1, (2), k , M, for each k = 1, ..., M and the current reference n> N by the method of matrix beams (IMF), after which the desired signal parameters are found by the formulas
Figure 00000040
а далее определяют текущую разность фаз и частоту мод колебаний, необходимую для определения параметров потока, протекающего через расходомер, обычно массового расхода и плотности жидкости, по формуламand then determine the current phase difference and the frequency of the vibration modes, necessary to determine the parameters of the flow flowing through the flow meter, usually the mass flow rate and the density of the liquid, according to the formulas
Figure 00000041
,
Figure 00000041
,
при этом на этом шаге вычисление параметров мод колебаний для отсчета n заканчивается, и они передаются в систему обработки данных потока, после чего в указанную систему поступает очередной, n+1-й блок отсчетов, что эквивалентно сдвигу окна оценки на один такт дискретизации, и цикл измерений повторяется снова.in this case, at this step, the calculation of the parameters of the oscillation modes for the n-count ends, and they are transferred to the flow data processing system, after which the next, n + 1-th block of samples arrives at the specified system, which is equivalent to shifting the estimation window by one sampling clock, and the measurement cycle repeats again. 2. Способ оценки текущей разности фаз и частоты сигналов y1(t) и у2(1) от 1-го и 2-го датчиков расходомера, представляемых суммой комплексных экспонент2. A method for estimating the current phase difference and frequency of signals y 1 (t) and 2 (1) from the 1st and 2nd flowmeter sensors, represented by the sum of complex exponents
Figure 00000042
Figure 00000042
где нижний цифровой индекс обозначает номер датчика, R1(2),k - комплексная амплитуда, k-й гармоники соответствующего датчика, а ωk, αk - частота и затухание, определяемые полюсом zk k-й гармоники, одинаковым для обоих датчиков, М - число экспоненциальных компонент в исследуемом сигнале, ε1(2),k(t) - аддитивная шумовая компонента соответствующего сигнала, отличающийся тем, что полюса сигналов принимают идентичными z1k=z2k=zk, k=1…М, при этом для расчета полюсов сигнала используются объединенные матрицы, получение в соответствии с векторным методом матричных пучков, что приводит к получению для каждого текущего отсчета n>N идентичных значений полюсов zk,n, k=1…М, которые используются для расчета соответствующих комплексных амплитуд R1(2),k,n, после чего параметры сигнала определяются по формуламwhere the lower digital index denotes the sensor number, R 1 (2), k is the complex amplitude of the kth harmonic of the corresponding sensor, and ω k , α k is the frequency and attenuation determined by the pole z k of the kth harmonic, the same for both sensors , M is the number of exponential components in the signal under investigation, ε 1 (2), k (t) is the additive noise component of the corresponding signal, characterized in that the signal poles are taken to be identical z 1k = z 2k = z k , k = 1 ... M, in this case, to calculate the signal poles, the combined matrices are used, obtaining in accordance with the vector method ternary beams, which leads to obtaining for each current reference n> N identical pole values z k, n , k = 1 ... M, which are used to calculate the corresponding complex amplitudes R 1 (2), k, n , after which the signal parameters are determined according to the formulas
Figure 00000043
Figure 00000043
и далее значения ωk,n и ϕk,n1,k,n2,k,n записываются в систему обработки данных оценки параметров потока, после чего в систему поступает очередной, n+1-й блок отсчетов, что эквивалентно сдвигу окна оценки на один такт дискретизации, и цикл вычислений повторяется снова.and further, the values of ω k, n and ϕ k, n = ϕ 1, k, n2, k, n are written to the data processing system for estimating the flow parameters, after which the next, n + 1-th block of samples which is equivalent to shifting the evaluation window by one sampling cycle, and the calculation cycle is repeated again.
RU2019113186A 2019-04-26 2019-04-26 Method for calculating current difference of phase and frequency of signals of inertial flow meters (versions) RU2707576C1 (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
RU2019113186A RU2707576C1 (en) 2019-04-26 2019-04-26 Method for calculating current difference of phase and frequency of signals of inertial flow meters (versions)

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
RU2019113186A RU2707576C1 (en) 2019-04-26 2019-04-26 Method for calculating current difference of phase and frequency of signals of inertial flow meters (versions)

Related Parent Applications (1)

Application Number Title Priority Date Filing Date
RU2017146968A Division RU2687803C1 (en) 2017-12-28 2017-12-28 Method for calculating current phase difference and frequency of signals of coriolis flowmeters

Publications (1)

Publication Number Publication Date
RU2707576C1 true RU2707576C1 (en) 2019-11-28

Family

ID=68836262

Family Applications (1)

Application Number Title Priority Date Filing Date
RU2019113186A RU2707576C1 (en) 2019-04-26 2019-04-26 Method for calculating current difference of phase and frequency of signals of inertial flow meters (versions)

Country Status (1)

Country Link
RU (1) RU2707576C1 (en)

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
RU2762219C1 (en) * 2021-01-11 2021-12-16 Олег Валентинович Жиляев Method for measuring the phase shift of signals from a coriolis flow meter
CN114235072A (en) * 2021-12-17 2022-03-25 电子科技大学 Zero-crossing detection-based Coriolis flowmeter phase difference calculation method

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US5578764A (en) * 1994-07-11 1996-11-26 Yokogawa Electric Corporation Coriolis type mass flowmeter utilizing phase shifters for phase shifting of the output signals
CA2593089A1 (en) * 2004-12-29 2006-07-06 Micro Motion, Inc. High speed frequency and phase estimation for flow meters
RU2448330C1 (en) * 2009-02-06 2012-04-20 Овал Корпорейшн Signal processing method, signal processing apparatus and coriolis acceleration flow metre
RU2460974C2 (en) * 2009-02-06 2012-09-10 Овал Корпорейшн Signal processing method, signal processing apparatus and coriolis flow meter

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US5578764A (en) * 1994-07-11 1996-11-26 Yokogawa Electric Corporation Coriolis type mass flowmeter utilizing phase shifters for phase shifting of the output signals
CA2593089A1 (en) * 2004-12-29 2006-07-06 Micro Motion, Inc. High speed frequency and phase estimation for flow meters
RU2371678C2 (en) * 2004-12-29 2009-10-27 Майкро Моушн, Инк. High-speed evaluation of frequency and phase of flow metres
RU2448330C1 (en) * 2009-02-06 2012-04-20 Овал Корпорейшн Signal processing method, signal processing apparatus and coriolis acceleration flow metre
RU2460974C2 (en) * 2009-02-06 2012-09-10 Овал Корпорейшн Signal processing method, signal processing apparatus and coriolis flow meter

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
RU2762219C1 (en) * 2021-01-11 2021-12-16 Олег Валентинович Жиляев Method for measuring the phase shift of signals from a coriolis flow meter
CN114235072A (en) * 2021-12-17 2022-03-25 电子科技大学 Zero-crossing detection-based Coriolis flowmeter phase difference calculation method
CN114235072B (en) * 2021-12-17 2023-04-18 电子科技大学 Zero-crossing detection-based Coriolis flowmeter phase difference calculation method

Similar Documents

Publication Publication Date Title
US7865318B2 (en) Meter electronics and methods for verification diagnostics for a flow meter
EP1949045B1 (en) Meter electronics and methods for determining one or more of a stiffness coefficient or a mass coefficient
US20100011882A1 (en) Method for operating a vibratory measuring instrument, and corresponding instrument
RU2707576C1 (en) Method for calculating current difference of phase and frequency of signals of inertial flow meters (versions)
US20130338943A1 (en) Method for operating a resonance measuring system and a resonance measuring system in this regard
Takamoto et al. New measurement method for very low liquid flow rates using ultrasound
Zhu et al. Mathematical modeling of ultrasonic gas flow meter based on experimental data in three steps
Henry et al. Response of a Coriolis mass flow meter to step changes in flow rate
CN101819056A (en) Instrument electronic device for checking and diagnosing flow meter and method
Thomsen et al. Perturbation-based prediction of vibration phase shift along fluid-conveying pipes due to Coriolis forces, nonuniformity, and nonlinearity
De Klerk et al. Uncertainty propagation in experimental dynamic substructuring
RU2515129C1 (en) Vortex flow meter
RU2687803C1 (en) Method for calculating current phase difference and frequency of signals of coriolis flowmeters
Rensing et al. Coriolis flowmeter verification via embedded modal analysis
JP7170049B2 (en) Method and apparatus for monitoring dissolution
JP2006184258A (en) Ultrasonic method and ultrasonic apparatus for computing concentration
Jamróz Effect of the continuous traverse trajectory and dynamic error of the vane anemometer on the accuracy of average velocity measurements at the cross-section of the mine heading–model-based testing
JP2616730B2 (en) Sound source motion analyzer
Semenov et al. Novel Prony-based algorithm for estimating oscillation parameters of Coriolis flowmeter at two-phase flow
Bause et al. An improved mode-tracing algorithm to compute dispersion curves of acoustic waveguides
RU2803043C1 (en) Method for assessing the state of a coriolis flowmeter for its verification and/or disagnostics
RU2403578C2 (en) Method of average flow rate determination
RU2762219C1 (en) Method for measuring the phase shift of signals from a coriolis flow meter
RU2773633C1 (en) Method for assessing the state of the measuring system of a coriolis flow meter
US11740117B2 (en) Method of determining total prove time

Legal Events

Date Code Title Description
MM4A The patent is invalid due to non-payment of fees

Effective date: 20201229