CC BY
33
УДК 53.082.22
ОСНОВАННЫЙ НА ФИЛЬТРЕ КАЛМАНА МЕТОД ИЗМЕРЕНИЯ ПЛОТНОСТЕЙ, ГРАНИЦЫ РАЗДЕЛА МЕЖДУ СЛОЯМИ И УРОВНЯ ДВУХСЛОЙНОЙ ЖИДКОСТИ В РЕЗЕРВУАРЕ
Е.Н. Абросимов, А.Л. Шестаков
BASED ON KALMAN FILTERING METHOD OF MEASUREMENT OF THE DENSITIES, THE INTERFACE BETWEEN LAYERS AND THE LEVEL OF TWO-LAYER LIQUID IN THE TANK
E.N. Abrosimov, A.L. Shestakov
Предлагается метод измерения плотностей, границы раздела менаду слоями и уровня двухслойной жидкости, основанный на применении фильтра Калмана. Приводятся результаты моделирования работы предлагаемого алгоритма.
Ключевые слова: измерение плотности, измерение уровня, определение положения границы между слоями, фильтр Калмана.
The article proposes a method of measuring densities, interface position between layers and level of two-layer liquid, based on application of Kalman filtering. Simulation results of the proposed algorithm are given.
Keywords: density measurement, level measurement, position estimation of interface between layers, the Kalman filter.
Введение
В ряде практических приложений (отстой товарной воды в нефтехранилищах, пенная сепарация и другие) возникает задача измерения уровня, положения границы раздела между слоями и плотностей двухслойных жидкостей.
В работе [1] был предложен метод измерения плотностей слоев, границы раздела между слоями и уровня двухслойной жидкости в резервуаре, главным достоинством которого является комплексное измерение требуемых параметров двухслойной жидкости. Одним из недостатков предложенного метода является отсутствие фильтрации шумов, которые неизбежно присутствуют в данных с применяемых датчиков давления.
Поэтому для повышения точности работы метода измерения предлагается использовать рекуррентный алгоритм, основанный на фильтре Калмана [2].
1. Начальные условия
Рассмотрим резервуар, в котором находится двухслойная жидкость, плотности р, и р2 слоев которой различаются существенно, то есть имеется ярко выраженная граница раздела между слоя-
Абросимов Евгений Николаевич ассистент кафедры информационно-измерительной техники, Южно-Уральский государственный университет аЬгозтюуИНп^.зши.ас.ги
Шестаков Александр Леонидович д-р техн. наук, профессор, ректор, Южно-Уральский государственный университет аёгшпй'игс. ас .ги
ми. В резервуар опущена штанга, на которой на расстоянии 2(1 закреплены два датчика давления 01 и £>2, находящиеся в разных слоях жидкости (рис. 1). Будем считать, что штанга движется по заданному закону х(ї). Примем, что начало координат находится в верхней граничной плоскости резервуара.
Рис. 1. Резервуар с датчиками шления на штанге
В работе [1] показано, что плотности слоев, положение границы раздела и общий уровень жидкости определяются следующим образом:
Abrosimov Evgeniy Nikolaevich - assistant lecturer of the Equipment for information and measuring department of South Ural State University; abrosimov@init.susu.ac.ru
Shestakov Aleksandr Leonidovich - PhD, professor, rector of South Ural State University; admin@urc.ac.ru
Р1(*2)~Р1(*1) .
д[х^2)-х^ 1)]’
Р2^2)~Р2^1) . д[х{Ь2)-х{Ь 1)]’
Р^-Ры
Нтр = х{12) + ар-^
Р13
р2а2)-р1а2)
(1) (2) (3)
, (4)
Р2-Р1 Я(Р2“Р1)
где Рл и Р2 - данные верхнего и нижнего датчика соответственно; Нж - общий уровень жидкости; Нгр - положение границы раздела между слоями;
- момент времени, когда нижний датчик давления находится у границы раздела, но не пересек ее; £2 - момент времени, когда верхний датчик давления находится у границы раздела, но не пересек ее.
При измерении в данных с датчиков давления будет неизменно присутствовать шум, который будет оказывать негативное влияние на точность измерения параметров рх, р2, Нж и Н, р. В связи с этим возникает задача повышения точности измерения по данным с датчиков давления О, и 02. Эту задачу можно успешно решить применением рекуррентного алгоритма, основанного на фильтрации Калмана.
Формулы (3) и (4) имеют явный нелинейный характер, что существенно усложняет применение фильтра Калмана. Для линеаризации задачи фильтрации произведем замену пары переменных Нж и Нгр на новую пару г = рхЯж и <7 = Я^Срг - рх), тогда формулы (3) и (4) будут выглядеть следующим образом:
^2) Ю (£2) Рп м)-
(5)
(6)
Определим связи меяеду переменными. Предварительно обозначим х(£2) — х(£,) = Дх. Из выражений (1) и (2) получим:
Р^2) - РЛк) =дЬх -рх; (7)
Р2&2)-Р2(11)=дАх-р2. (8)
Из выражения (5) можно записать:
А выражение (6) перепишем следующим образом:
Полученные уравнения запишем в систему:
+д(х(Ь) + сО ■ р2 + 0 ■ г + (-д) ■ ц, которая в матричном виде запишется следующим образом:
(11)
~ ^(А) Р2(£2) — ^(а) Р\ (^г) — Р.п м Р2(£2) — (^г)-
дАх О
(12)
Эта матричная запись справедлива в условиях отсутствия шумов.
2. Стохастические уравнения, описывающие ход технологического процесса
Для использования фильтра Калмана необходимо описать технологический процесс в виде разностных уравнений, которые дают представление о нормальном ходе процесса. Отклонение от нормального хода процесса принимается случайным и считается шумом процесса.
Будем считать, что типичное поведение параметров ръ р2, г и (7 - медленное изменение, приближенно линейное на коротком интервале времени. Используя подход, предложенный в работе [3], запишем уравнения для параметров в разностном виде: ^
откуда
Р1,к+1 = 2р 1'к ~ Р1,к-1 шР1,(с- (13а)
Аналогично
откуда
Р2./С + 1 = 2р2,к - Р2,к-1 + ШР2Х (136)
Для параметра г также можно записать:
откуда
2/С + 1 — 2Zk ~ %к-1 "I" (13в)
Аналогичным образом для параметра <7 получим:
откуда
<7/с+1 = 2Як ~ + шЧк. (13г)
где к - индекс, нумерующий дискретные моменты
времени; сор | /;. оор2 /;. оо2/; и шумы процесса.
Примем, что шумы процесса имеют нормальное распределение с нулевыми математическими ожи-
7 7 7 7
Даниями И дисперсиями <5р1к, ®р2к> и °с[кС00Т~ ветственно.
Запишем уравнения (13а)-(13г) в матричном виде. Для этого составим систему:
(14)
2к’
Запишем ее в матричном виде:
РіЛ+і'
Рі ,к ?2,к + 1 ?2 ,к гк +1
Чк + 1 Чк
2-10 0 О О О О
10000000 002 -1 0000 00100000 00002 -1 00 00001000 0000002 -1 Ц) 0000010-
Рі.к ] ГШР1*
РіЛ-1 о ?2,к Р2,/с-1 ?к 2к-1 Як
-І I- О
X
+
<А
Р2,к
О
од.
"гк О
(15)
Введем обозначения:
-Рі,/с+1 Рі ,к Рг,/с+1 Р2 ,к гк+1 2/с
Чк + 1 ■ Чк
времени /і'+1:
- вектор состояния в момент
1 О
О О
О О
- мат-
00001000
0000002 -1 Ц) 0 0 0 0 0 10-
рица эволюции вектора состояния Хк в отсутствии шумов;
РіЛ Р 2,к
времени к;
Рг,/с—і ?к 2к-1 Чк
<7/с-1
вектор состояния в момент
Од
Рі./С
о
Р2./С
о
Од
о
Од
- вектор шума процесса.
Чк
О
В результате система (14) запишется следующим образом:
Хк+1=А-Хк + УУк. (16)
3. Уравнения, описывающие процесс измерения
Уравнение (16) - это уравнение процесса. Добавим к нему уравнение измерения, полученное из выражения (12) с использованием введенного вектора Хк и с учетом шумов измерения:
Рі,к(*2) ~ р1,к&і) = 3&х ' Рі,к + 0 ■ 9і,к-і
+д Ах ■
Р2,к-1 + 0 ■ гк + 0 ■ гк_! О ■ о^_1 + 1 ■ гг.
РцеС^) Ра™ —
(17)
Г2,к($2) Рі,к(і2) —
Чк-і + 1 ■ у2к\
О ■ аь_, + 1 ■ 17„
+ ° ' Р1,к-1 + 0 '
+0 ■ гк_1
р2,
В матричном виде система (16) может быть записана следующим образом:
^1,/с(^2) — РглИг)
^г./сС^г) — Рг.к^О _
Рг,к(р2) ~ Рцтш
?2,к(р2^ ~ Р1,кИ2)-
дкх о ° о о о о о'
о о д&с о о о о о
д(х(£2)-Ю о о о -д о о о
-д(х&2) - Ю о д(х(с2) + ^) о о о -д о.
хХк +
(18)
Введем следующие обозначения:
Рі,к0-2) ~ Рі,к(Рі)
Рі.кО-г) ~ ^атм
РгкІЇї) Рі .кО-г)-
вектор результатов
измерения;
о
д(х(і2) - (і)
С =
О
дАх
О
О О О О
о -д
-д(х(с2) -сі) о д(х&2) ■ матрица наблюдения;
а) о
о -д о.
- вектор шума измерения.
В итоге система (17) в матричном виде запишется следующим образом:
Ук = С-Хк + Ук. (19)
В результате получили динамическую стохастическую систему:
(20)
4. Синтез фильтра Калмана
На первом этапе синтеза фильтра задаются начальные значения параметров вектора состояния Х0 и матрицы ковариации Р0. Вектор состояния в начальный момент времени:
Р1,0"|
Р1,0 Р2,0 Р2,0
20 Чо Чо
Этот вектор состоит из величин рх 0, р2 0. г0 и <70, вычисленных по данным датчиков давления, полученным на начальном интервале движения:
^1 (^2,о) ^1 (^1,о) .
(21)
1,0 Я[ж(‘2,0)-ж(‘і,0)]’
_ Р2(^2,о)~Р2,о(^1,о). '° я[ж(‘2,о)-ж(‘і,о)].
(22)
(23)
(24)
20 — Рі,0^ж,(Ь Чо = ^гр,о(Рг,о — Рі,о), (25)
где А о и і2,о _ моменты времени на начальном
интервале движения, когда нижнии и верхнии (соответственно) датчики давления находятся у границы раздела, но не пересекли ее.
Примем, что в начальный момент времени штанга была погружена в жидкость так, что датчики О, и 02 находились в разных слоях жидкости.
Определимся с начальными значениями элементов ковариационной матрицы ошибок оценивания вектора состояния в начальный момент вре-
мени. В общем виде ковариационная матрица Р0 определяется выражением [2]:
Р0 = м[(Х0-Х0)(Х0-Х0)Т1 (26)
где Х0 - истинное значение вектора состояния в начальный момент времени. Элементами вектора Х0 — Х0 будут являться среднеквадратические отклонения погрешностей определения плотностей слоев в начальный момент времени. Тогда ковариационная матрица Р0 запишется в виде
где а,
Рі,0
0
0
о
о
о
о
о
2
Р 1,0’
о
°о
Р 1,0
о о о о о о 12
~>Р 2,0’
о
о
ао
Р2.0
о
о
о
о
о
т2
о
о
о
°о
Р2.0
о
о
о
о
т2
о
о
о
о
сг|
20
о
о
о
о
о
о
о
о
о
<
о
о
о
о
о
о
о
о
<-1
(27)
°г„ и 0дГ1 априорные значения
дисперсий в начальный момент времени. При этом \ 2
'Чо
, (р2,0 Рі,о)
т2 _ _2 „2 >г0 ~ оНГРі0Рі,0
■ нтр/(о:
грЛ
т2
2
Р 1,0
Р 2,0
(28)
С ростом числа итераций к матрица ковариации вектора состояния Рк быстро стабилизируется и перестает зависеть от своего начального значения.
Задав начальные значения, можно сделать прогноз ожидаемых значений матриц фильтра Калмана. Прогнозируем вектор состояния Хк по Хк_1 и ковариационную матрицу Рк по Рк-±:
Хк+і =А-Хк; (29)
Рк+1 — А
Рк-АТ + (1,
(30)
где - ковариационная матрица шума процесса, которая характеризует влияние шумов друг на друга, то есть характеризует отклонение процесса от линейности. Эта матрица определяется следующим образом:
'Рі,к
О
о
о
о
о
о
о
о
о
*Рг,к
о
о
о
о
о
о
о
о
о
<
о
о
о
Так как шумы процесса сор
О
О
О
О
О
О
О
О
(л)„
о
о
о
о
о
о
о
(31)
Р1,к’ Р2,к
независимы, все недиагональные элементы матрицы <2 равны нулю.
Следующим этапом синтеза фильтра является корректировка прогнозируемых значений. Для этого необходимо вычислить оптимальную матрицу усиления фильтра (коэффициент Калмана) [2]:
ик=Рк-Ст -[С-Рк-Ст + Д]-\ (32)
где Я - ковариационная матрица шума измерения, которая характеризует влияние погрешностей датчиков давления друг на друга. Матрица Л определится следующим образом:
^k 0 ^ О 2 ст? О
где
а\
О
(ЗЗ)
>1 ^ 2J
дисперсии погрешностей датчиков давления. Считаем что погрешности датчиков давления не коррелированы между собой.
Обновим оценку для среднего значения вектора состояния по новым результатам измерения Ук:
Хк=Хк + ик-[¥к-С-Хк], (34)
и оценку для матрицы ковариации:
Рк = Рк-ик-С-Рк. (35)
В процессе работы фильтра с увеличением числа итераций к матрицы Рк и ик стабилизируются, что говорит о сходимости фильтра Калмана. Стоит заметить, что предельные значения этих матриц будут зависеть от матриц Л, О и А.
5. Численное моделирование работы
алгоритма и фильтра Калмана
Предложенный алгоритм был смоделирован в пакете МаНаЬ 7.7.0 со следующими параметрами:
- время моделирования 20 ООО с;
- высота резервуара 5 м;
- плотность верхнего слоя 150 кг/м3 и за время моделирования увеличивается до 350 кг/м3;
- плотность нижнего слоя 800 кг/м3 и за время моделирования увеличивается до 1200 кг/м3;
- в начальный момент времени граница раздела находилась на расстоянии 3,4м от верхней граничной плоскости резервуара, в течение первых 10 000 с уровень границы раздела снизился до 3,6 м, а затем за вторые 10 000 с повысился до 3,4 м;
- в начальный момент времени уровень двухслойной жидкости находился на расстоянии 1,8 м от верхней граничной плоскости резервуара, в течение первых 10 000 с уровень повысился до 1,2 м, а затем за вторые 10 000 с снизился до 1,8 м;
- предполагалось, что погрешности датчиков представляют собой белый шум с амплитудой 2 %
от предельного значения показания датчиков с учетом атмосферного давления.
Для простоты моделирования начальное значение координаты штанги было принято равным координате положения границы раздела между слоями.
На рис. 2-5 представлены графики результатов моделирования.
По результатам моделирования определили, что:
- СКО плотности р, верхнего слоя
api= 4,73кг/м3;
- СКО плотности р, верхнего слоя с использованием фильтра Калмана
api =0,21 кг/м3;
- СКО плотности р2 нижнего слоя
ар 2= 14,44кг/м3;
- СКО плотности р2 нижнего слоя с использованием фильтра Калмана
Р2
= 0,87 кг/м3;
- СКО положения границы раздела
ан = 3,78-10 м;
пгр ? ’
- СКО положения границы раздела с использованием фильтра Калмана
1 м;
аЪ = 1,38-10"4
нгп ?
- СКО уровня двухслойной жидкости ан = 1.54 • 10 3 м;
- СКО уровня двухслойной жидкости с использованием фильтра Калмана
of, = 1,87-10 м.
Мж ?
Значения СКО an
'Р1- ~Р2- и V, Рассчитывались по простому методу, предложенному в [1]. Для всех параметров смещение математического ожидания отсутствует.
Значения СКО и сТр2 брались из диагональных элементов с координатами [ 1; 1] и [3;3]
матрицы кова-
соответственно установившейся риации Р. Значения СКО ст^ и а^ж рассчитывались по формулам (13) и (15).
кис. 2. Плотность pt верхнего слоя
О 50 100 150 200 250 300 350 400
Нис. 3. Плотность р2 нижнего слоя
Рис. 4. Положение границы раздела НТ„ между слоями
Рис. 5. Уровень нж двухслойной жидкости
Результаты проведенного моделирования показывают, что применение рекуррентного алгоритма, основанного на фильтре Калмана, приводит к существенному повышению точности измерения.
Заключение
В статье предложен метод измерения плотностей, положения границы раздела между слоями и уровня двухслойной жидкости в резервуаре, основанный на применении фильтрации Калмана. Приведены результаты моделирования предложенного метода, которые подтверждают его работоспособность и более высокую точность по сравнению с расчетами этих параметров прямым методом.
Литература
1. Абросимов, Е.Н. Измерение плотностей, границы раздела между слоями и общего уровня двухслойной жидкости в резервуаре / Е.Н. Абросимов, А.Л. Шестаков // Вестник ЮУрГУ. Серия «Компьютерные технологии, управление, радиоэлектроника». -2011. -Вып. 14. -№23(240). - С. 4-8.
2. Браммер, К. Фильтр Калмана-Бьюси. Детерминированное наблюдение и стохастическая фильтрация / К. Браммер, Г. Зиффлинг. - М.: Наука, 1982.
3. Применение фильтра Калмана в задаче измерения уровня и плотности жидкости с помощью двух датчиков давления / О.Л. Пбряева, П.Г. Ко-репанов, А.С. Семенов, А.Л. Шестаков // Информационно-измерительная техника. - 2007. - № 6.
Поступила в редакцию 7 апреля 2011 г.
