Способ вычисления текущей разности фаз и частоты сигналов кориолисовых расходомеров (варианты)
METHOD FOR CALCULATING CURRENT DIFFERENCE OF PHASE AND FREQUENCY OF SIGNALS OF INERTIAL FLOW METERS (VERSIONS)

Изобретение относится к контрольно-измерительной технике и способам обработки одного или более сигналов датчиков в расходомере и может быть использовано в приборостроении при разработке и изготовлении кориолисовых расходомеров.
Массовые расходомеры Кориолиса широко используются для измерения массового расхода плотности и объемного расхода, а также получения другой информации о веществах, протекающим через трубопровод, как раскрыто в патенте US 4491025, МПК G01F 1/76, G01F 1/84, опубл. 01.01.1985 г. Эти расходомеры имеют одну или более расходомерных труб различных конфигураций. Каждая конфигурация трубы имеет набор мод собственных колебаний, в том числе, но не только, плоский изгиб, крутильную, радиальную и связанную моду. В типичном варианте применения измерения массового расхода методом Кориолиса используемая конфигурация труб возбуждается на одной или более колебательных мод при протекании вещества через трубопровод, и движение расходомерных труб измеряется в точках, разнесенных по длине трубы, как показано на фиг. 1
Приводной механизм возбуждает колебания расходомерной трубки. Когда нет вещества, протекающего через расходомер, все точки вдоль расходомерной трубки колеблются с идентичной фазой. По мере того как вещество начинает протекать через расходомерную трубку, ускорения Кориолиса приводят к тому, что каждая точка вдоль расходомерной трубки имеет различную фазу относительно других точек вдоль расходомерной трубки. Фаза на входной стороне расходомерной трубки запаздывает от приводного механизма, тогда как фаза на выходной стороне опережает приводной механизм.
Датчики помещаются в различных точках на расходомерной трубке и преобразуют движение этих точек в соответствующий набор сигналов. Разность фаз соответствующих мод сигналов датчиков сигналов пропорциональна массовому расходу вещества, протекающего через расходомерную трубку или расходомерные трубки.
Из предыдущего уровня техники известны способы обработки сигналов для измерения массового расхода Кориолиса, основные из которых рассмотрены в [1]. Для оценки разности фаз необходимой моды колебаний используются дискретное преобразование Фурье (DFT), быстрое преобразование Фурье (FFT), различные модификации режекторных и следящих фильтров. Далее полученные оценки разности фаз и частоты колебаний блока расходомерных трубок могут быть использованы для того, чтобы вычислить массовый расход и плотность расходуемого вещества.
Указанные способы работают достаточно хорошо в стационарном режиме, когда расходуемое вещество в расходомере является однородным и значения мод колебаний, определяемые плотностью и расходом контролируемого вещества, стабильны.
Однако, когда расходуемое вещество является неоднородным, например, в двухфазных потоках, где расходуемое вещество содержит жидкость и твердое вещество или имеются пузырьки газа в жидком расходуемом веществе, возникают быстрые флюктуации частот мод колебаний, которые не могут отслеживаться техникой предыдущего уровня, в силу принципиальных ограничений преобразования Фурье и основанных на нем фильтров.
Это эффект также является проблемой в режимах «пустой - полный - пустой», когда возможно многократное появление значительных объемов газа эквивалентных временному полному исчезновению измеряемой среды из расходомерных трубок.
Описанные режимы приводят к некорректным измерениям частоты, что, в свою очередь влечет значительные погрешности (20% и более) в измерении разности фаз, а значит, в дальнейшем, плотности и расхода.
Наиболее близким к заявляемому решению является способ, предложенный для решения указанной проблемы в патенте RU 2371678, МПК G01F 1/84 опубл. 27.10.2009 «Высокоскоростная оценка частоты и фазы расходомеров» и заключающийся в использовании преобразования Гильберта, реализованном на основе 90 - градусных фазовращателях. Способ обработки сигналов датчиков в расходомере для последующего вычисления массового расхода плотности и объемного расхода реализуется электронным измерительным оборудованием, содержащим интерфейс для приема первого сигнала датчика и второго сигнала датчика и связанную с интерфейсом систему обработки, предназначенную для формирования девяностоградусного сдвига фаз из первого сигнала датчика с помощью преобразования Гилберта и вычисления разности фаз из девяностоградусного сдвига фаз, первого сигнала датчика и второго сигнала датчика. Частоту вычисляют из первого сигнала датчика и девяностоградусного сдвига фаз. Второй девяностоградусный сдвиг фаз может быть сформирован из второго сигнала датчика.
Способ позволяет существенно повысить динамические характеристики расходомера и уменьшить погрешность оценки разности фаз в условиях небольших и достаточно медленных вариации частоты.
Тем не менее, при существенных и быстрых изменениях частоты колебаний расходомерных трубок, что характерно для режимов «полный - пустой - полный» или скачков плотности расходуемого вещества (резкое изменение газовой фазы в двухфазных средах), вышеперечисленные недостатки методов предшествующего уровня обнаруживает и указанный способ. Это связано с тем, что реализация широкополосных фазовращателей базируется на тех же методах обработки сигналов, что и методы предыдущего уровня техники, а значит и имеет те же принципиальные ограничения по динамике измерений, что и они. Кроме того, преобразование Гильберта ориентировано на работу с одной модой колебаний, что не выполняется при многофазном потоке. Для уменьшения влияния шумовых компонент в приведенном способе предлагается использовать фильтрацию, что дополнительно ухудшает быстродействие алгоритма и снижает эффективность его использования.
Несмотря на различные способы оценки разности фаз и частоты, приводящие к улучшению отдельных характеристик кориолисовых расходомеров, работающих с многофазными средами, все они, в том числе и использованный в качестве прототипа, имеют один и тот же недостаток. Эти преобразования можно отнести к непараметрическим методам [2], имеющими принципиальное ограничение на разрешение частот, связанное с временем наблюдения примерным соотношением Δω≈1/ΔT, где Δω и ΔT - необходимое разрешения по частоте и время наблюдения необходимое для его обеспечения, соответственно. Это соотношение накладывает жесткие требования на длительность наблюдаемого участка при требованиях повышенного разрешения, что в свою очередь ухудшает динамические характеристики алгоритмов обработки и затрудняет работу с нестационарными сигналами.
Задача, на решение которой в первую очередь направлено заявляемое изобретение, заключается в уменьшении погрешности измерения массового и объемного расхода жидкой среды при наличии возмущающей фазы (газовой или твердой) в кориолисовых расходомерах за счет сокращения количества отсчетов и, соответственно, времени, требуемого для получения текущего значения расхода.
Указанная задача достигается тем, что согласно первому варианту сигналы от датчиков положения расходомерных трубок кориолисова расходомера представляются в виде суммы комплексных экспонент
1(2),k
1(2),k
1(2),k
1(2),k
где нижний цифровой индекс обозначает номер датчика, R, ω, α - комплексная амплитуда, частота и затухание k-ой гармоники соответствующего датчика,  - полюса сигналов, М - число экспоненциальных компонент в исследуемом сигнале, ε(t) - аддитивная шумовая компонента соответствующего сигнала. Для оценки неизвестных параметров до начала вычислений, на основе известного диапазона частот колебаний и динамических характеристик измерительной системы, определяют исходные параметры: период дискретизации сигналов Т, число отсчетов N, определяющее длительность текущего окна обработки и некоторый целочисленный параметр L,  и сохраняют их в системе обработки данных. Входные сигналы с 1-го и 2-го датчиков расходомера подают на указанную систему обработки, где сигналы дискретизируют с заданным периодом Т, и преобразуют в дискретную последовательность отсчетов от 1-го и 2-го датчиков
1(2),k
1(2),k
k
k
n=0…S>N, n - номер текущего отсчета, которые используют для определения параметров сигнала первого и второго датчиков R, z, М, для каждого k=1, …, М и текущего отсчета n>N. Для вычисления параметров R, z, М используют метод матричных пучков (ММП) [3], после чего вычисляют искомые параметры сигнала по формулам
Далее определяют текущую разность фаз и частоту мод колебаний, необходимую для определения параметров потока, протекающего через расходомер, обычно массового расхода и плотности жидкости, соответствующие номеру n текущего отсчета. При этом частота моды и разность фаз по первому способу определяется по формулам (4)
После этого, окно оценки параметра сдвигается на один дискрет, т.е. n1=n+1 и в качестве нового окна оценки используется новый массив отсчетов n1-N÷n1=n+1-N÷n+1 и вышеприведенная последовательность действий повторяется.
1
2
Согласно второму варианту, сигнал y(t) и y(t) от 1-го и 2-го датчиков расходомера, представляется суммой комплексных экспонент
1(2),k
k
k
k
1(2),k
1k
2k
k
k,n
1(2),k,n
где нижний цифровой индекс обозначает номер датчика, R - комплексная амплитуда, k-ой гармоники соответствующего датчика, а ω, α - частота и затухание, определяемые полюсом z k-ой гармоники, одинаковым для обоих датчиков, М - число экспоненциальных компонент в исследуемом сигнале, ε(t) - аддитивная шумовая компонента соответствующего сигнала, согласно изобретению изначально предполагается идентичность полюсов сигналов, т.е. z=z=z, k=1..М, при этом, для расчета полюсов сигнала используются объединенные матрицы, получение в соответствии с векторным методом матричных пучков, что приводит к получению для каждого текущего отсчета n>N идентичных значений полюсов z, k=1..М, которые используйся для расчета соответствующих комплексных амплитуд R, после чего параметры сигнала определяются по формулам
k,n
k,n
1,k,n
2,k,n
и далее значения ω и ϕ=ϕ-ϕ записываются в систему обработки данных оценки параметров потока, после чего в систему поступает очередной, n+1 - й блок отсчетов, что эквивалентно сдвигу окна оценки на один такт дискретизации, и цикл вычислений повторяется снова.
Уменьшение времени обработки происходит за счет того, что в параметрических методах минимальное число отсчетов определяется не требуемым разрешением по частоте, как в непараметрических методах, а числом неизвестных параметров, которое может быть в десятки и сотни раз меньше, чем число отсчетов на интервале разрешения частоты методами Фурье и другими подобными методам. При этом, в первом приближении, длительность дискреты в параметрическом методе не регламентируется, т.е. формально необходимое число отсчетов может быть получено за интервал меньший, чем период колебании самой высокочастотной моды, в то время как для непараметрических методов необходимый интервал наблюдения может составлять десятки периодов. Учитывая это, предлагаемые способы выдают некоторую «точечную» оценку параметров в окне (на интервале наблюдения), получение зависимости этих сигналов от времени реализуется путем смещения окна (интервала наблюдения) при получении нового отсчета. Алгоритм смещения окна при получении нового отсчета отображают соответствующие шаги, показанные на фиг. 2 и 3.
Сущность изобретения поясняется следующими графическими материалами:
Фиг. 1 Блок-схема кориолисова расходомера, где позициями обозначены следующие элементы: корпус расходомера 1, расходомерные трубки 2, датчики 3;
Фиг. 2 Блок-схема алгоритма по 1 варианту способа;
Фиг. 3 Блок-схема алгоритма по 2 варианту способа;
Фиг. 4 Сравнение алгоритмов в одномодовом режиме работы; отображены результаты расчетов по первому варианту способа и с использованием преобразования Гильберта по патенту RU 2371678, а также фактическое значение измеряемых параметров. Для демонстрации зависимости от длительности окна оценки графики и расчеты приведены для N=25 и N=100. В соответствующей таблице приведены оценки погрешности методов при измерении частоты, разности фаз и амплитуды сигналов.
Фиг. 5 Сравнение алгоритмов при наличии второй моды; отображены результаты расчетов по второму варианту способа и с использованием преобразования Гильберта по патенту RU 2371678, а также фактическое значение измеряемых параметров при наличии дополнительной синусоидальной компоненты, являющейся в данном случае помехой, т.е моделируется двухмодовый режим работы. Добавленная синусоида имеет частоту 100 Гц. И амплитуду 20 dB от основной. Расчеты погрешностей методов для оцениваемых параметров приведены в таблице.
Фиг. 6 Оценки и истинное значения полюсов мод, рассчитанные при помощи варианта 1 (верхняя и средняя диаграммы) и варианта 2 (нижняя диаграмма). Истинное значение полюсов z_1=0,8070+0,5863j и z_2=0,8484+0,5199j отмечено кружочком, а оцененное крестиком. Концентрация оценок вокруг истинных значений на нижних рисунках свидетельствуете об уменьшении погрешности при использовании второго варианта способа, что подтверждают данные приведенные в соответствующей таблице.
Фиг. 7 Частота сигналов, рассчитанная по варианту 1 (сверху и посередине) и варианту 2 (снизу). Отметим значительное уменьшение дисперсии оценки частоты на нижнем графике (вариант 2);
Фиг. 8 Амплитуда, рассчитанная по варианту 1 (светлым) и варианту 2 (темным), для правой и левой катушки (сигнал 1, сигнал 2). Видно значительное уменьшение дисперсии оценки амплитуды при использовании второго варианта способа;
Фиг. 9 Разница фаз сигналов, вычисленная по способу 1 (синим) и способу 2 (красным). Отметим уменьшение дисперсии оценки разности фаз, что свидетельствует об уменьшении погрешности при использовании способа 2.
В качестве основы модели сигнала расходомера, работающего в условиях двухфазного потока использовалась модель, предложенная в [2].
При двухмодовом режиме работы мощность «мешающей» моды принималось 20 dB от основной.
Алгоритм заявляемого способа по первому варианту приведен на фиг. 2 и содержит следующие шаги:
1. Установка исходных параметров алгоритма, их сохранение в вычислительной системе, дискретизация сигналов и передача каждого отсчета по мере его получения в вычислительную систему.
2. Формирование массива («окна») из N>=4M последовательных отсчетов n-М÷n по каждому датчику, где n - текущий отсчет.
a
b
3. Формирование матриц Y, Y (5), (6) для каждого датчика.
a
4. Нахождение SVD-разложения матриц Y.
5. Оценивание число М полюсов сигнала.
6. Нахождение усеченных до ранга М псевдообратных матриц  (10).
k
E
7. Оценивание полюсов z для каждого датчика путем вычисления собственных значений матриц Z (12).
k
8. Оценивание амплитуд R для каждого датчика, путем решения матричного уравнения (13) методом наименьших квадратов.
9. Вычисление параметров мод в соответствии с (14), (15).
10. Смещение «окна» на одну позицию, т.е. n=n+1.
На первом шаге способа устанавливаются необходимый период дискретизации сигналов Т и размер «окна» обработки N>4M. После их сохранения в вычислительной системе с помощью АЦП реализуется процесс дискретизации и получения отсчетов сигнала, которые по мере их поступления также передаются в вычислительную систему.
k
k
Далее, в соответствии с ММП, находят полюсы z с помощью решения задачи на обобщенные собственные значения пучка матриц, составленных из отсчетов сигнала у(n). По найденным полюсам ММП находят оценки для R и параметры мод колебаний. Изложим подробнее способ с использованием ММП.
a
b
После дискретизации сигналов с датчиков и формирования массивов измерений по каждому датчику (шаг 2), формируют две матрицы Y, Y (шаг 3) размеров (N-L)×L, следующим образом:
k
Здесь М≤L≤N-М параметр метода. Показано [3], что при выборе  дисперсия оценки полюсов z будет минимальна, т.е. ММП будет наименее чувствителен к шуму. Для конкретизации можно положить  (функция Round округляет число до ближайшего целого).
Основой метода матричных пучков является сингулярное разложение (SVD разложение).
a
b
Для матриц Y, Y справедлива следующая факторизация [3]:
Здесь
b
a
L
R
k
b
a
k
Рассмотрим пучок матриц Y-λY=ZR(Z-λE)Z, где Е - единичная матрица -го порядка. Можно показать [3], что при М≤L≤N-М и λ≠z, k=1, 2, …, М, ранг матрицы Y-λY равен М. Однако, если λ=z, то k-ая строка матрицы Z-λЕ - нулевая и ранг этой матрицы равен М-1.
b
a
Таким образом, полюсы  могут быть найдены как обобщенные собственные значения пучка матриц Y-λY [3], т.е. как собственные значения матрицы  (Здесь индекс  используется для обозначения псевдообратной матрицы.)
a
В случае зашумленных данных матрицу Y для эффективной фильтрации шума необходимо подвергнуть операции сингулярного (SVD) разложения
a
Здесь U, V - унитарные матрицы, S - диагональная матрица, элементами которой являются сингулярные числа матрицы Y. (Верхний индекс Т означает операцию транспонирования.)
Заметим, что в случае отсутствия шума диагональная матрица S имеет ровно М ненулевых сингулярных чисел, все последующие равны нулю. В случае зашумленного сигнала ненулевых сингулярных чисел уже не будет, однако между первыми М и последующими сингулярными числами матрицы S будет наблюдаться ярко выраженный скачок, который и позволит определить число комплексных экспонент в сигнале.
a
Сингулярное разложение матрицы Y позволяет определить число истинных экспонент сигнала М. Кроме того, оно может быть использовано для нахождения псевдообратной матрицы  На практике используется усеченная до ранга М псевдообратная матрица:
1
M
0
1
M
где σ, …, σ - это М наибольших сингулярных чисел матрицы  соответствующие им сингулярные векторы, , S=diag(σ, …, σ).
k
После нахождения матрицы  для оценки полюсов z, k=1, …, М, сигнала остается только найти М собственных чисел матрицы  или, в силу следующей цепочки равенств:
найти М собственных чисел матрицы
k
k
Далее, в методе матричных пучков по известным М и z находятся комплексные амплитуды R из решения следующей задачи:
k
После вычисления комплексных амплитуд R можно определить искомые параметры сигнала по формулам
Далее, на шаге 9 вычисляют текущую разность фаз и частоту мод колебаний, необходимую для определения параметров потока, протекающего через расходомер, обычно массового расхода и плотности жидкости, по формулам
На этом шаге вычисление параметров мод колебаний для отсчета n заканчивается, и они передаются в систему оценки параметров потока, после чего на шаге 10 в систему поступает очередной, n+1 - й блок отсчетов, что эквивалентно сдвигу окна оценки на один такт дискретизации.
Модификации метода и иные способы реализации ММП не изменяют сущности варианта.
Согласно второму варианту способ реализуется путем использования информации об идентичности полюсов сигнала с датчиков, предложенным в [5].
Алгоритм второго варианта способа приведен на фиг. 3 и содержит следующие шаги:
1. Установка исходных параметров алгоритма, их сохранение в вычислительной системе, дискретизация сигналов и передача каждого отсчета по мере его получения в вычислительную систему.
2. Формирование массива («окна») из N>=4M последовательных отсчетов n-М÷n по каждому датчику, где n-текущий отсчет.
a
1
b1
a
2
b2
3. Формирование матриц Y, Y Y, Y (5), (6) для каждого датчика.
4. Формирование объединенных матриц пучка
a
U
5. Нахождение SVD-разложения матрицы Y.
6. Оценивание числа М полюсов сигнала.
7. Нахождение усеченной до ранга М объединенной псевдообратной матрицы  (10).
k
E
8. Оценивание М объединенных полюсов z путем вычисления собственных значений матрицы Z, полученной на основе (12) с использованием усеченной матрицы
k
k
9. Оценивание амплитуд R для каждого датчика, путем решения матричного уравнения (13) методом наименьших квадратов с подстановкой значений полюсов z, полученных на предыдущем шаге.
10. Вычисление параметров мод в соответствии с (14), (15) (с учетом получения параметров частоты и затухания по объединенной матрице).
11. Смещение «окна» на одну позицию, т.е. n=n+1.
Основным отличием второго варианта от первого варианта способа является вычисление общих полюсов сигнала по объединенным матрицам пучка. В связи с этим первые три шага варианта совпадают с первыми тремя шагами основного способа.
k
bU
a
U
a
b
a
U
bU
mk
Далее, в целях получения оценок объединенных полюсов составляются блочные матрицы пучка  Как показано в [5], М полюсов z являются обобщенными собственными значениями пучка матриц Y-λY. В связи с этим значения объединенных полюсов могут быть получены в соответствии с алгоритмами первого варианта способа путем замены матриц Y, Y на матрицы Y, Y. То есть шаги 5, 6, 7, 8 второго варианта аналогичны шагам 4, 5, 6, 7 первого варианта способа при соответствующей замене матриц. Комплексные амплитуды R на шаге 9 алгоритма находятся затем из решения 2 систем уравнений:
1
M
ik
k
где z…z - объединенные полюса, полученные на шаге 8 алгоритма. Вычисление интересующих параметров мод по полученным значениям R и z, i=1..M, k=1, 2 проводится на шаге 10 по выражениям (14), (15) аналогично первому варианту способа, с учетом идентичности полюсов для первого и второго сигнала.
Для сравнения характеристик предлагаемого алгоритма и прототипа было проведено математическое моделирование их работы. Для генерации исходных сигналов была использована модель, предложенная в [4] со следующими параметрами: основная частота изменяется в пределах 85-100 Гц, амплитуда основной частоты изменяется в пределах 0.05-0.3 В, разность фаз - в пределах 0-4 град., частота дискретизации 2 кГц. Для устранения резких изменений параметров их реализации были пропущены через низкочастный фильтр с полосой 6 Гц.
Выводы. Из приведенных результатов расчетов следует, что предложенный способ решает задачу по уменьшению погрешности измерения частоты и разности фаз при малых временах измерения, причем как показывает моделирование при уменьшении времени измерения преимущество предлагаемого метода увеличивается. Наличие дополнительных мод сокращает разницу между методами, но и в этом случае преимущество предлагаемого метода остается существенным. Уменьшение погрешности измерения частоты и разности фаз позволяют предположить, что соответствующим образом будет уменьшена и погрешность оценки расхода и плотности контролируемой жидкости при двухфазных режимах работы. Учет информации об идентичности полюсов согласно варианта 2 может существенно уменьшить погрешность измерения параметров. Однако при этом необходимо учитывать возможное усложнение системы обработки за счет увеличения размерности матриц, а также конкретные варианты конструкции расходомера, возможно приводящие к не идентичности потоков в расходомерных трубках.
Список литературы
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.
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.
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.