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 PDFInfo
- 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
Links
Images
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01F—MEASURING VOLUME, VOLUME FLOW, MASS FLOW OR LIQUID LEVEL; METERING BY VOLUME
- G01F1/00—Measuring the volume flow or mass flow of fluid or fluent solid material wherein the fluid passes through a meter in a continuous flow
- G01F1/76—Devices for measuring mass flow of a fluid or a fluent solid material
- G01F1/78—Direct mass flowmeters
- G01F1/80—Direct mass flowmeters operating by measuring pressure, force, momentum, or frequency of a fluid flow to which a rotational movement has been imparted
- G01F1/84—Coriolis or gyroscopic mass flowmeters
Abstract
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
Приводной механизм возбуждает колебания расходомерной трубки. Когда нет вещества, протекающего через расходомер, все точки вдоль расходомерной трубки колеблются с идентичной фазой. По мере того как вещество начинает протекать через расходомерную трубку, ускорения Кориолиса приводят к тому, что каждая точка вдоль расходомерной трубки имеет различную фазу относительно других точек вдоль расходомерной трубки. Фаза на входной стороне расходомерной трубки запаздывает от приводного механизма, тогда как фаза на выходной стороне опережает приводной механизм.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
Способ позволяет существенно повысить динамические характеристики расходомера и уменьшить погрешность оценки разности фаз в условиях небольших и достаточно медленных вариации частоты.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
где нижний цифровой индекс обозначает номер датчика, R1(2),k, ω1(2),k, α1(2),k - комплексная амплитуда, частота и затухание k-ой гармоники соответствующего датчика, - полюса сигналов, М - число экспоненциальных компонент в исследуемом сигнале, ε1(2),k(t) - аддитивная шумовая компонента соответствующего сигнала. Для оценки неизвестных параметров до начала вычислений, на основе известного диапазона частот колебаний и динамических характеристик измерительной системы, определяют исходные параметры: период дискретизации сигналов Т, число отсчетов N, определяющее длительность текущего окна обработки и некоторый целочисленный параметр L, и сохраняют их в системе обработки данных. Входные сигналы с 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, 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, 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
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
Далее определяют текущую разность фаз и частоту мод колебаний, необходимую для определения параметров потока, протекающего через расходомер, обычно массового расхода и плотности жидкости, соответствующие номеру 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)
После этого, окно оценки параметра сдвигается на один дискрет, т.е. 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
где нижний цифровой индекс обозначает номер датчика, 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
и далее значения ωk,n и ϕk,n=ϕ1,k,n-ϕ2,k,n записываются в систему обработки данных оценки параметров потока, после чего в систему поступает очередной, n+1 - й блок отсчетов, что эквивалентно сдвигу окна оценки на один такт дискретизации, и цикл вычислений повторяется снова.and further, the values of ω k, n and ϕ k, n = ϕ 1, k, n -ϕ 2, 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:
Фиг. 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 (
Фиг. 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
В качестве основы модели сигнала расходомера, работающего в условиях двухфазного потока использовалась модель, предложенная в [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. Нахождение усеченных до ранга М псевдообратных матриц (10).6. Finding truncated to rank M pseudoinverse matrices (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:
Здесь М≤L≤N-М параметр метода. Показано [3], что при выборе дисперсия оценки полюсов zk будет минимальна, т.е. ММП будет наименее чувствителен к шуму. Для конкретизации можно положить (функция Round округляет число до ближайшего целого).Here M≤L≤N-M parameter of the method. It is shown [3] that when choosing 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 (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]:
ЗдесьHere
Рассмотрим пучок матриц 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.
Таким образом, полюсы могут быть найдены как обобщенные собственные значения пучка матриц Yb-λY a [3], т.е. как собственные значения матрицы (Здесь индекс используется для обозначения псевдообратной матрицы.)So the poles can be found as generalized eigenvalues of the matrix pencil Y b -λY a [3], ie as eigenvalues of a matrix (Here is the index 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
Здесь 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 позволяет определить число истинных экспонент сигнала М. Кроме того, оно может быть использовано для нахождения псевдообратной матрицы На практике используется усеченная до ранга М псевдообратная матрица: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 In practice, a pseudoinverse matrix truncated to rank M is used:
где σ1, …, σM - это М наибольших сингулярных чисел матрицы соответствующие им сингулярные векторы, , S0=diag(σ1, …, σM).where σ 1 , ..., σ M is the M largest singular numbers of the matrix their corresponding singular vectors, , S 0 = diag (σ 1 , ..., σ M ).
После нахождения матрицы для оценки полюсов zk, k=1, …, М, сигнала остается только найти М собственных чисел матрицы или, в силу следующей цепочки равенств:After finding the matrix to estimate the poles z k , k = 1, ..., M, the signal can only be found M eigenvalues of the matrix or, by virtue of the following chain of equalities:
найти М собственных чисел матрицыfind M eigenvalues of the matrix
Далее, в методе матричных пучков по известным М и 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:
После вычисления комплексных амплитуд Rk можно определить искомые параметры сигнала по формуламAfter calculating the complex amplitudes R k, it is possible to determine the desired signal parameters by the formulas
Далее, на шаге 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
На этом шаге вычисление параметров мод колебаний для отсчета 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. Формирование объединенных матриц пучка 4. Formation of combined beam matrices
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. Нахождение усеченной до ранга М объединенной псевдообратной матрицы (10).7. Finding a truncated to rank M joint pseudoinverse matrix (10).
8. Оценивание М объединенных полюсов zk путем вычисления собственных значений матрицы ZE, полученной на основе (12) с использованием усеченной матрицы 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
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.
Далее, в целях получения оценок объединенных полюсов составляются блочные матрицы пучка Как показано в [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 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
где 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
Для сравнения характеристик предлагаемого алгоритма и прототипа было проведено математическое моделирование их работы. Для генерации исходных сигналов была использована модель, предложенная в [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.,
Выводы. Из приведенных результатов расчетов следует, что предложенный способ решает задачу по уменьшению погрешности измерения частоты и разности фаз при малых временах измерения, причем как показывает моделирование при уменьшении времени измерения преимущество предлагаемого метода увеличивается. Наличие дополнительных мод сокращает разницу между методами, но и в этом случае преимущество предлагаемого метода остается существенным. Уменьшение погрешности измерения частоты и разности фаз позволяют предположить, что соответствующим образом будет уменьшена и погрешность оценки расхода и плотности контролируемой жидкости при двухфазных режимах работы. Учет информации об идентичности полюсов согласно варианта 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
Список литературы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)
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)
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)
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 |
-
2019
- 2019-04-26 RU RU2019113186A patent/RU2707576C1/en not_active IP Right Cessation
Patent Citations (5)
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)
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 |