CC BY
8
Вычислительные технологии
Том 22, № 5, 2017
Исследование частотно-избирательных свойств преобразования Гильберта — Хуанга и его модификаций на примере изучения автоколебательных пульсаций давления
А. А. Левин, В. А. Спиряев*
Институт систем энергетики им. Л. А. Мелентьева СО РАН, Иркутск, Россия *Контактный e-mail: errolorr@gmail.com
Рассмотрено применение преобразования Гильберта — Хуанга (ПГХ) к анализу временных рядов, порожденных автоколебательными пульсациями давления в объеме вскипающей жидкости. Сформулировано понятие качественной декомпозиции, полученной с помощью ПГХ. Показано принципиальное отличие с точки зрения качества декомпозиции базового ПГХ от его модификации с применением адаптивного шума. Дополнительно проанализированы различные параметры модифицированного ПГХ.
Ключевые слова: преобразование Гильберта —Хуанга, автоколебательные пульсации давления, нестационарный временной ряд.
Введение
В настоящее время существует необходимость обеспечения устойчивой к возможным возмущениям и безопасной работы теплообменных устройств различного назначения в энергетической, химической и микроэлектронной технике, характеризующейся высокими удельными тепловыми нагрузками. Актуальность экспериментальных методов исследования протекающих сложных тепловых, гидродинамических процессов в условиях возникновения кризиса теплоотдачи обусловливается многообразием возможных сочетаний теплофизических параметров и геометрических факторов. В работах [1-3] представлены результаты экспериментального изучения процессов образования и эволюции паровой структуры на теплоотдающих тонкопроволочных и пленочных нагревателях.
Естественным ограничителем для жидкостных теплоносителей является наступление кризиса теплоотдачи, зачастую протекающего в условиях резко нестационарного роста температуры теплоотдающей поверхности. Один из важных аспектов изучения развития нестационарных кризисных состояний — определение динамики давления в объеме вскипающей жидкости. Межфазные превращения, протекающие в условиях охлаждения поверхности, перегретой выше температуры насыщения жидкости, зачастую сопровождаются возникновением импульсов давления различной интенсивности.
© ИВТ СО РАН, 2017
В случае реализации термоакустического эффекта, впервые описанного Рэлеем, возможна генерация автоколебательных пульсаций давления.
В работе [4] выявлены особенности генерации автоколебаний при кипении недогре-тых жидкостей в трубах. Отмечается малая изученность автоколебательных процессов, протекающих при кипении жидкости в узких кольцевых каналах в условиях наступления кризиса теплоотдачи. В работе [5] авторы указывают на совпадение частот пульсаций центров парообразования с частотой пульсации парожидкостной среды как необходимое условие для возникновения резонансных термоакустических автоколебаний, создающих стоячие волны большой амплитуды. Особое внимание автоколебательным процессам в двухфазных потоках уделено в обзоре [6].
На практике изучение автоколебательных пульсаций давления сводится к анализу временных рядов, ими порожденных. Такие временные ряды описывают нестационарные сигналы, амплитудные и частотные характеристики которых очень быстро изменяются во времени. Поэтому для повышения качества их анализа необходимо использовать подход, обладающий свойством адаптивности к каждому конкретному сигналу. Одним из эффективных методов для исследования и анализа таких нестационарных временных рядов является преобразование Гильберта — Хуанга (ПГХ).
В работе [7] Хуангом предложен отличный от традиционного подход к декомпозиции временных рядов, основанный не на априорном выборе базисных функций, а на их адаптивном построении в процессе обработки входной информации. В отличие от классических методов, для которых признаком останова является получение остатка в виде белого шума, построение базисных функций в методе Хуанга, напротив, начинается с самой высокочастотной составляющей, а признаком останова является остаток в виде монотонной функции. На втором этапе осуществляется построение аналитической функции, мнимая компонента которой связана с полученной аппроксимацией исходного вещественного временного ряда интегральным преобразованием Гильберта, позволяющим получить мгновенные частоты и амплитуды конкретных составляющих. Этим объясняется закрепившееся в литературе название ПГХ. Таким образом, можно говорить, что полученные составляющие совместно со своими частотами и амплитудами составляют модель исходного процесса.
Изложенная выше простая и эффективная идея ПГХ находит применение в огромном числе различных по своей природе прикладных задач. Так, например, в [7] рассмотрено применение ПГХ для анализа волн, производимых ветром в океане. Другим применением ПГХ является использование его для идентификации различных режимов двухфазного перехода газа в твердое состояние [8]. В работе [9] ПГХ используется в качестве метода предобработки данных для опознавания источников беспроводных излучений. Представляют интерес возможности ПГХ в задачах выделения тренда и фильтрации шума для нестационарных сигналов. Например, в [10] исследуются измеренные с помощью PMU (Phasor Measurement Units) колебания большой объединенной энергетической системы для выявления и устранения шумовых составляющих и трендов. В работах [11-14] ПГХ применяется в задачах геофизики, медико-биологических исследованиях, в энергетике, обработке изображений и речевых сигналов. Необходимо отметить работы французского коллектива авторов [15-18], посвященные исследованию теории и практическим возможностям ПГХ.
В данной статье рассматривается вопрос использования ПГХ и его модификации для адаптивной обработки временных рядов, описывающих динамику давления в объеме вскипающей жидкости, и анализа частотно-избирательных свойств этих методов.
1. Алгоритмы ПГХ и модифицированного ПГХ
Рассмотрим алгоритмы классического ПГХ (далее просто ПГХ) и его модифицированной версии (мПГХ).
ПГХ состоит из двух независимых частей: метода разложения на эмпирические моды (Empirical Mode Decomposition — EMD) и преобразования Гильберта. Основная идея EMD состоит в разложении исследуемого сигнала x(t) на эмпирические базисные функции с последующим применением к ним преобразования Гильберта. Название таких функций в русскоязычной литературе еще не устоялось. Их называют эмпирическими модами, модальными функциями, характеристическими функциями, базисными функциями и даже, основываясь на дословном переводе английского названия, внутренними модовыми функциями. Все эти названия представляются не совсем точными и удачными, поэтому в дальнейшем будем использовать сокращение от английского названия — IMF (Intrinsic Mode Function). Процедура декомпозиции x(t), вообще говоря, представляет собой итеративный процесс, который удобнее всего представить в виде блок-схемы (рис. 1).
Как видно из рис. 1, алгоритм EMD подразумевает два цикла: внутренний цикл "отсеивания" и внешний цикл декомпозиции. На внутреннем цикле выделяются конкретные IMF, а на внешнем определяется их количество и производится завершение работы алгоритма, если выполнен выбранный критерий останова. Каждая полученная IMF должна удовлетворять двум следующим свойствам:
• число экстремумов и число пересечений оси абсцисс на заданном интервале должны совпадать;
• для любого t среднее значение, определяемое верхней и нижней огибающими, должно быть близко к нулю.
В данной работе критерий останова алгоритма на внутреннем уровне обеспечивается вышеуказанными свойствами IMF, хотя существуют и другие критерии останова [7,17,20]. Останов на внешнем уровне осуществляется, когда количество экстремумов в последней IMF меньше двух. Такая IMF называется остатком разложения.
После работы алгоритма EMD любой сигнал x(t) может быть представлен в следующем виде:
где q < р < n, п — общее количество выделенных IMF, из которых c\(t),. . . , cq (t) — высокочастотные компоненты, cq+\(t),... ,ср(t) — компоненты, определяющие физические свойства ряда, и cp+\(t),... , cn-\(t), rn(t) — трендовые и несинусоидальные компоненты, при этом rn(t) является остатком разложения.
Метод EMD, входящий в ПГХ, помимо своих очевидных достоинств, таких как адаптивность и зависимость только от исходных данных, имеет и недостатки, связанные с наличием в одной IMF осцилляций, слишком различных по значениям амплитуды, или распределением по нескольким IMF одной осцилляции. Оба этих недостатка получили название "смешивание мод". Для решения данной проблемы предложен метод множественной декомпозиции — Ensemble EMD (EEMD) [21], основанный на многократном добавлении к сигналу белого шума заданной амплитуды и вычислении средних значений всех полученных IMF для выполнения итогового разложения. Согласно
(1)
Конец
Рис. 1. Блок-схема метода разложения на IMF
свойствам EMD [16], предложенная процедура позволяет решить проблему "смешивания мод", однако имеет свои недостатки, основные из которых — неодинаковое количество IMF в разложениях с различным шумом и большая вычислительная сложность. Предложенная в [22, 23] вариация EEMD позволяет, сохраняя преимущества EEMD, преодолеть и его недостатки.
Кратко рассмотрим эту вариацию, получившую название CEEMDAN (Complete Ensemble Empirical Mode Decomposition with Adaptiv Noise). Определим оператор Ej(■), который для заданного сигнала производит j-ю моду с помощью EMD. Пусть шг есть белый шум с N(0,1). Тогда для заданного сигнала x(t) метод CEEMDAN можно описать с помощью следующего алгоритма:
1) применим обычный EMD к /-реализациям x(t)+ e0^z(t), г = 1,I, и получим первую моду
— 1 1
i м f (t) = jJ2IMFl (t);
i=1
2) для к = 1 подсчитаем первый остаток
n(t) = x(t) - IMF1(t);
3) разложим реализации r1(t) + е1Е1(шг(t)), г = 1,I, до их первой моды и определим вторую
— 1 1
I M F2(t) = -^Е^ n(t)+ е1Е1(шг(Ш
г=1
4) для к = 2,..., К вычислим к-й остаток
rk(t) = rk-1(t) - IMFk(t);
5) так же, как и в шаге 3, найдем ( к + 1)-ю моду
— 1 1
IMFk+1(t) = jY, Е1( rk (t) + екЕк (ujl(t)));
г=1
6) повторяем шаг 4 для следующих к, пока не получим остаток, который удовлетворяет заданному критерию останова. В этом случае алгоритм останавливается, а последний остаток имеет вид
к
R(t) = x(t) - ^ 1ММк(t), к=1
где К — количество полученных мод. Таким образом, исходный сигнал x(t) можно представить в следующем виде:
к
x(t) = ^ TMFk (t) + R(t). (2)
k=1
Формула (2), в отличие от EEMD, делает предложенную декомпозицию завершенной и дает точную реконструкцию исходного сигнала.
Отметим, что коэффициенты е¿, вообще говоря, позволяют на каждом этапе выбирать разную амплитуду добавляемого белого шума, измеряемую в терминах амплитуды исходного временного ряда или соотношения сигнал/шум в децибелах. Однако в данной работе такая возможность не рассматривалась. В работе [21] предложено использовать
маленькие амплитуды шума для сигналов, находящихся преимущественно в высокочастотной области, и наоборот.
Следующий этап в преобразовании Гильберта — Хуанга состоит в применении к каждой IMF преобразования Гильберта. Это позволяет выделить из данной IMF две составляющие: амплитудную и фазовую. Рассмотрим это более подробно для произвольного действительного сигнала x(t).
Дополним сигнал x(t) до аналитической функции
z(t) = x{t) + ixH(t), (3)
где г — мнимая единица; Xu(t) — преобразование Гильберта, определяемое как
xn(t) = 1 P.V. / ^ds. (4)
к J t - s
— <х
Интеграл в (4) понимается как главное (P.V. — principal value) значение по Коши. Перепишем (3) в экспоненциальной форме
z (t) = A(t)etm, (5)
где
A(t) = у! x2(t) + xH(t)
мгновенная амплитуда, а
т = и<*ё (7)
— фаза исходного сигнала х(Ь). Тогда, согласно определению, производная от (7) определяет мгновенную частоту
ф) = ф(1) = | агоЬЕ ^. (8)
С учетом (6) и (8) исходный сигнал можно представить в следующей форме:
if uij (s)ds I
£Aj(t)e 0 j.
, _ . if шл (s)ds
x(t) = ReaW > A3(t)e 0 }. (9)
В (9) остаток rn может быть отброшен, а может быть включен в зависимости от задачи и физических предположений. Отметим, что в представление (9) амплитуда и частота входят как функции времени и называются мгновенными амплитудой и частотой. В силу этого факта представление сигнала через набор IMF принципиально отличается от представления Фурье
х(Ь) = И,еа1 < ^^ а^е
где аз и являются константами. Таким образом, ПГХ представляет собой обобщенное преобразование Фурье. Формула (9) позволяет представить мгновенные амплитуды
и частоты как функции времени в трехмерном пространстве. Это частотно-временное распределение амплитуд определяется как спектр Гильберта
п 4
i f Шу (s)ds
=1
H(u, t) = ^Aj(t)e 0
Дополнительно к спектру Гильберта Хуангом [7] определен маргинальный спектр
т
к(ш) = ! И(и, 5 )(1з, (10)
0
который представляет собой общий вклад амплитуды, соответствующей каждой мгновенной частоте на протяжении всей длины сигнала, и является аналогом фурье-спектра.
2. Экспериментальные расчеты
При изучении временных рядов теплофизической природы с помощью ПГХ проблема "смешивания мод" довольно сильно проявляется и существенно затрудняет анализ полученных составляющих. Использование мПГХ позволяет в той или иной степени решить эту проблему.
Рассмотрим один ряд из серии временных рядов, которые представляют собой автоколебания давления, возникающие в процессе импульсного тепловыделения в недогре-той до температуры насыщения жидкости. Чередующиеся процессы вскипания жидкости и последующей конденсации переохлажденного водяного пара порождают стоячие волны в канале с собственной частотой, близкой к 48 Гц. По мере нестационарного охлаждения стенки интенсивность циклических фазовых переходов снижается, что приводит к изменению амплитудно-частотных характеристик пульсаций. Одним из важных факторов для понимания механизмов взаимодействия жидкости, пара и перегретой стенки в этих условиях является масштаб мелкочастотных пульсаций межфазной поверхности, свидетельствующих о возможности реализаций режима высокоинтенсивного теплообмена. Наиболее доступным способом анализа динамики межфазной поверхности является измерение давления в жидкости с применением высокочастотных датчиков. Полученные ряды данных, один из них представлен на рис. 2, в дальнейшем обрабатывались методами ПГХ и мПГХ. Анализ динамики межфазной поверхности заключается в изучении полученных после декомпозиции характеристических функций (в дальнейшем IMF) и их частот.
Исследование проводилось в пакете прикладных программ для решения технических и математических задач MATLAB R2013b с использованием алгоритмов ПГХ и мПГХ [24].
Согласно методу мПГХ, для разложения исходного сигнала необходимо задать три параметра: амплитуду белого шума ( N), размер ансамбля ^N) и количество итераций для работы классического метода EMD. Количество итераций для метода EMD не играет существенной роли и по сути служит удобным аналогом критерия останова метода на внутреннем уровне. В то же время N и Е N — важные параметры для получения качественной декомпозиции. Вообще говоря, оценку качества декомпозиции и важности
t
Рис. 2. Автоколебание давления
полученных IMF можно проводить по-разному. Например, в [25] предлагается использовать статистические оценки, позволяющие установить, насколько энергия конкретного IMF больше энергии белого шума. Такие оценки позволяют только определить, является ли конкретная IMF белым шумом или нет, но совершенно не касаются проблемы качества самих IMF. Поэтому в работе под качественной декомпозицией будем понимать такое разложение исследуемого сигнала, которое, по возможности, состоит из "несмешанных" IMF. При этом будем предполагать, что каждая "несмешанная" IMF соответствует одному из физических процессов. Соответствие конкретной IMF тому или иному процессу определяется особенностями самого исследуемого сигнала и экспертом. Демонстрацию качественной и некачественной декомпозиции с помощью метода EMD и семейства методов CEEMDAN с различными параметрами N и EN покажем на уже выбранном сигнале (рис. 2).
В результате разложения сигнала методами EMD и CEEMDAN с различными значениями N и EN получены различные наборы IMF, содержащие 12-15 составляющих и монотонный остаток. Оценка качества набора заключается в анализе частотного состава IMF, входящих в этот набор. Для определения "смешанности" IMF будем использовать метод оценки спектральной плотности мощности (СПМ). Существует несколько алгоритмов для оценки СПМ сигнала, основанных на использовании только значений исходного временного ряда: метод периодограмм, методы модифицированных периодограмм (Даньелла, Бартлетта, Уэлча) и метод Блэкмана — Тьюки. В настоящей работе используется периодограмма Уэлча, так как ее оценка СПМ является асимптотически несмещенной и состоятельной и, как следствие, менее осциллирующей, чем у других методов.
Для начала покажем принципиальное отличие методов EMD и CEEMDAN с разумно выбранными параметрами амплитуды шума и размера ансамбля, например, как предлагается в [22], N = 0.02 и EN = 500. После декомпозиции с помощью предложенных методов получим два набора IMF, представленных на рис. 3.
Для удобства разделим оба набора на три группы по частоте с одинаковым количеством составляющих для обоих методов. Первая группа состоит из высокочастотных составляющих (IMF1-IMF5), вторая группа — из средне- и низкочастотных (IMF6-IMF8), а низкочастотная группа (IMF9-IMF14) имеет составляющие с очень низкими амплитудами, которые не несут важной информации о сигнале и могут быть отброшены. В эту группу попадает и вторая IMF, полученная с помощью CEEMDAN, которая тоже имеет слишком низкую амплитуду и не влияет на формирование исходного сигнала. Приме-
0.05 0
-0.05
600 700 800 900 1000 1100 1200
-н-н-
0 0 кллха кШкккШ i IL L L L L LLLL LlfJ
600 700 800 900 1000 1100 1200
0.05 F 0 •
-0.05 t
'НИМИ
0.05 0
-0.05
52 0.05 -
e 0" -0.05 L
-H4-hmw
.ffffffHf.
55 °.°5 S „ „0
0.005 с 0 ■
-0.005 t
0.04 С 0 -
-0.041
0.04 f
1 4-.H ri/i -
ю 0.02 f
¡3 0 -
к -0
4 H>">....... *...........M.....HHf frfif"
нн>-
t»»»
I 0
" -0.2
S 0
6 0.2 1 0 Д -0.2
900 1000 1100 t, мс
900 1000 t, мс
1100 1200
700
800
600
700
800
1200
Рис. 3. Наборы IMF для методов EMD (слева) и CEEMDAN (справа)
нение периодограммы Уэльча к первым двум группам IMF позволяет сравнить СПМ двух методов.
Для большей наглядности покажем СПМ только для составляющих из второй группы. Из графиков на рис. 4 хорошо видно, что метод CEEMDAN позволяет разделить составляющие с достаточно явно выраженной частотой в отличие от EMD, составляющие которого имеют гораздо более размазанную структуру и совпадающую частоту. Для первой (высокочастотной) группы применение СПМ дает аналогичный результат. Отметим, что получение "несмешанного" (в той или иной степени) набора IMF возможно для достаточно широкого разброса параметров, влияющих на работу алгоритма CEEMDAN. Это показывает, что метод CEEMDAN позволяет существенно лучше проводить декомпозицию исходного сигнала для определения частот конкретных IMF. Однако, как видно из рис. 4, IMF все же смешиваются. Так, IMF7 имеет две ярко выраженных частоты, а IMF6 вообще размазана по всему частотному спектру, поэтому естественно подобрать наилучшие значения амплитуды шума и размера ансамбля для получения наиболее "несмешанного" набора составляющих (IMF). Для оценки качества IMF будем использовать метод оценки СПМ.
Рассмотрим варианты выбора параметров разложения метода CEEMDAN для того же сигнала. Согласно рекомендациям [22], EN достаточно выбрать равным 500. В то же время результаты расчетов, приведенные в [19], показывают, что уже при EN > 100 декомпозиция сигналов различной природы является вполне приемлемой (наборы практически не отличаются друг от друга) и не требует повышения размера ансамбля. Однако для сигналов, исследуемых в работе, выбор EN оказался чувствителен к амплитуде шума N: чем выше N, тем больше надо выбирать EN. Поэтому для качественного усреднения при выделении высокочастотных составляющих сигнала в дальнейших расчетах выбран параметр EN = 500. Дальнейшее увеличение размера ансамбля не дает заметного результата, но при этом требует гораздо больших вычислительных ресурсов. Более существенным в получении качественной декомпозиции является изменение параметра амплитуды белого шума N. Рассмотрим этот параметр более подробно.
Для сигнала, представленного на рис. 2, проведена декомпозиция при разных значениях N. Рассмотрены разложения с отклонением шума в диапазоне от 0.002 до 0.3 с шагом 0.002. Так как основная энергия сигнала содержится в шестой, седьмой и восьмой IMF, в дальнейшем будем рассматривать на предмет смешанности именно эти составляющие. Отметим, что для значений N < 0.018 декомпозиция напоминает декомпозицию,
2.5 2 1.5 1
0.5 0
PSD10
—IMF6 —IMF7 .......IMF8
4 4.5 5 5.5
log2(ffl), Гц
6.5 7 7.5
PSD10
2.521.5 1
0.5
-IMF6 -IMF7 IMF8
0
4 4.5 5 5.5
log2(ü>), Гц
6.5 7 7.5
Рис. 4. Спектральная плотность мощности группы среднечастотных IMF для методов EMD (a) и CEEMDAN (б)
а
полученную с помощью классического метода EMD. Это объясняется очень маленькими значениями N. При таких значениях происходит сильное смешивание одних IMF в другие, как показано на рис. 5, а. Уже с 0.018 < N < 0.06 IMF начинают разделяться (рис. 5, б), однако только со значения N = 0.06 IMF разделяются полностью. Максимально качественная декомпозиция получается при 0.066 < N < 0.08. В этом интервале происходит максимальное разделение и амплитуды разделенных IMF достигают максимальных значений. На рис. 5, б приведен СПМ для N = 0.07. При этом уже для N = 0.1 опять начинает происходить эффект "смешивания мод". На рис. 5, в хорошо видно, как энергия восьмой IMF частично перешла в седьмую, а сама седьмая IMF смешалась с шестой. Таким образом, при N = 0.07 получаем декомпозицию (рис. 5, б), в которой наиболее чисто выделяются составляющие с частотой около 44 и 85 Гц, соответствующие седьмой и восьмой IMF, и более зашумленная, составляющая с частотой около 114 Гц — шестая IMF. Отметим, что подбор амплитуды шума позволяет полностью очистить восьмую IMF и наиболее качественно и точно определить ее частоту, в то время как шестую и седьмую IMF не удалось полностью разделить, что может говорить о дефекте метода или действительном (физическом) содержании двух и более частот в этих составляющих. На рис. 6 представлены частоты шестой, седьмой и восьмой IMF.
а б в
2.5 2 1.5 1
0.5 0
PSD10
4 4.5 5 5.5
—IMF6 —IMF7 .......IMF8
.5 7 7.5
log2(ffl), Гц
2.5 2 1.5 1
0.5 0
PSD10
4 4.5 5 5.5
—IMF6 —IMF7 .......IMF8
log2(ffl), Гц
2.521.5 1
0.5
PSD10
.5 7 7.5
4 4.5 5 5.5
—IMF6 — IMF7 .......IMF8
log2(ffl), Гц
.5 7 7.5
Рис. 5. Спектральная плотность мощности для шестой, седьмой и восьмой IMF, полученных с помощью метода CEEMDAN при значениях амплитуды шума: N = 0.01 (a); N = 0.07 (б); N = 0.2 (в)
я
^ 120 ^ 80 $ 40 0
я
^ 120 S 80
£ 40
0
Я
^ 120-S 80£ 40; 0L
700
800
900 1000 t, мс
1100
1200
Рис. 6. Мгновенные частоты шестой, седьмой и восьмой IMF
£ 2000
^ 1600 §
900 950 1000 1050 1100
t, мс
Рис. 7. Мгновенная частота первой IMF
Подобный результат можно получить и для высокочастотных составляющих. И хотя для хорошего разделения первых IMF требуется практически отсутствие шума, результат не сильно искажается при EN > 500 и N = 0.07. Анализируя набор высокочастотных IMF, удалось выявить важную высокочастотную составляющую с частотой 1600-1700 Гц. Соответствующая ей первая IMF имеет сильно зашумленный вид, так как исходный сигнал содержит большой шумовой фон, поэтому выделить эту составляющую можно, только сопоставив ее с соответствующей ей амплитудой. На рис. 7 представлена частота первой IMF.
С физической точки зрения представленные на рис. 6 и 7 графики частот характеризуют механизмы автоколебательных пульсаций давления. С помощью предложенной методики удалось идентифицировать возникающие автоколебания со средней частотой 44 и 85 Гц, соответствующие седьмой и восьмой IMF. Анализ видеокадров из [26] также показал наличие стоячих волн давления с характерной скоростью 9 м/с в канале длиной 35 см, что соответствует частоте около 51 Гц. С усилением амплитуды пульсационно-го процесса в спектре идентифицируется высокочастотная составляющая колебаний давления с частотой 1600-1700 Гц, соответствующая первой IMF (рис. 7). Проявление этой составляющей совпадает с моментами достижения максимумов давления и позволяет идентифицировать две стадии автоколебательного процесса: начальную (с 450 до 800 мс) и стадию развитого кипения (с 850 до 1400 мс). Наступление стадии развитого кипения хорошо отслеживается и на низкочастотных IMF (см. рис. 6). Более подробно протекающие процессы и их сопоставление со скоростной видеосъемкой рассмотрены в работе [26].
Заключение
Рассмотрены метод ПГХ и его модификация для изучения автоколебательных пульсаций давления. На реальном сигнале, принятом за тестовый, отработаны методы ПГХ и мПГХ. Результаты исследования показали абсолютное преимущество метода мПГХ. Дополнительно получены параметры (амплитуда шума и размер ансамбля) мПГХ для проведения более качественной декомпозиции и вычисления соответствующих этой декомпозиции мгновенных частот и амплитуд.
Предложенная методика по выбору метода и его параметров может быть легко адаптируема для изучения и анализа других сигналов как такой же, так и иной природы, а при определенном совпадении амплитудных и частотных характеристик сигнала могут использоваться даже полученные в работе конкретные значения амплитуды шума и размера ансамбля. Использование правильных параметров при разложении исследуемого сигнала на составляющие позволяет более точно выделять компоненты, соответствующие тем или иным механизмам, порождающим исходный сигнал.
Благодарности. Работа выполнена на оборудовании уникальной научной установки ИСЭМ СО РАН "Высокотемпературный контур" при частичной финансовой поддержке РФФИ (грант № 15-01-01425).
Список литературы / References
[1] Glod, S., Poulikakos, D., Zhao, Z., Yadigarogly, G. An investigation of microscale explosive vaporization of water on an Ultrathin Pt Wire // Intern. J. of Heat Mass Transfer. 2002. Vol. 45. P. 367-379.
[2] Pavlenko, A.N., Chekhovich, V.Yu. Interconnection between dynamics of liquid Boiling-up and heat transfer crisis for nonstationary heat release // J. of Eng. Thermophys. 2007. Vol. 16, No. 3. P. 175-187.
[3] Surtaev, A.S., Pavlenko, A.N. Observation of boiling heat transfer and crisis phenomena in falling water film at transient heating // Intern. J. of Heat and Mass Transfer. 2014. Vol. 74, No. 7. P. 342-352.
[4] Дорофеев Б.М., Волкова В.И. Гидродинамические и термоакустические автоколебания при поверхностном кипении в каналах // Акустический журнал. 2008. Т. 54, № 5. С. 732-739.
Dorofeev, B.M., Volkova, V.I. Hydrodynamic and thermoacoustic self-oscillations at surface boiling in channels // Acoustical Phys. 2008. Vol. 54, No. 5. P. 633-639.
[5] Smirnov, H.F., Zrodnikov, V.V., Boshkoa, I.L. Thermoacoustic phenomena at boiling subcooled liquid in channels // Intern. J. Heat Mass Transfer. 1997. Vol. 40, No. 8. P. 19771983.
[6] Ruspini, L.C., Marcel, C.P., Clausse, A. Two-phase flow instabilities: A review // Intern. J. Heat Mass Transfer. 2014. Vol. 71, No. 4. P. 521-548.
[7] Huang, N.E., Shen, Z., Long, S.R., Wu, M.C. et al. The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis // Proc. of the Royal Soc. A. 1998. Vol. 454. P. 903-995.
[8] Hu, H.L., Zhang, J., Dong, J., Luo, Z.Y., Xu, T.M. Identification of gas-solid two-phase flow regimes using Hilbert — Huang transform and neural-network techniques // Instrum. Sci. & Technology. 2011. Vol. 39, No. 2. P. 198-210.
[9] Yuan, Y., Huang, Z., Wu, H., Wang, X. Specific emitter identification based on Hilbert — Huang transform-based time-frequency-energy distribution features // IET Communications. 2014. Vol. 8, No. 13. P. 2404-2412.
[10] Laila, D.S., Messina, A.R., Pal, B.C. A refined Hilbert — Huang transform with applications to interarea oscillation monitoring // IEEE Transact. on Power Sys. 2009. Vol. 24, No. 2. P. 610-620.
[11] Battista, B.M., Knapp, C., McGee, T., Goebel, V. Application of the empirical mode decomposition and Hilbert — Huang transform to seismic reflection data // Geophysics. 2007. Vol. 72, No. 2. P. H29-H37.
[12] Echeverria, J.C., Crowe, J.A., Woolfson, M.S., Hayes-Gill, B.R. Application of empirical mode decomposition to heart rate variability analysis // Medical and Biological Eng. and Comput. 2001. Vol. 39, No. 4. P. 471-479.
[13] Browne, T.J., Vittal, V., Heydt, G.T., Messina, A.R. A comparative assessment of two techniques for modal identification form power system measurements // IEEE Trans. on Power Sys. 2008. Vol. 23, No. 3. P. 1408-1415.
[14] Алимурадов А.К. Исследование частотно-избирательных свойств методов декомпозиции на эмпирические моды для оценки частоты основного тона речевых сигналов // Тр. Моск. физ.-техн. ин-та. 2015. Т. 7, № 3(27). С. 56-68.
Alimuradov, A.K. Research of frequency-selective properties of empirical mode decomposition methods for speech signals' pitch frequency estimation // Tr. Mosk. Fiz.-tekhn. In-ta. 2015. Vol. 7, No. 3(27). P. 56-68. (In Russ.)
[15] Flandrin, P., Goncalves, P., Rilling, G. Detrending and denoising with empirical mode decomposition // Proc. of the 12th Europ. Signal Proc. Conf. (EUSIPCO '04). Vienna: IEEE Publisher, 2004. Vol. 2. P. 1581-1584.
[16] Flandrin, P., Rilling, G., Goncalves, P. Empirical mode decomposition as a filter bank // IEEE Signal Proc. Letters. 2004. Vol. 11, No. 2. P. 112-114.
[17] Rilling, G., Flandrin, P., Go^alves, P. On empirical mode decomposition and its algorithms // IEEE-EURASIP Works. on Nonlinear Signal and Image Proc. NSIP-03, Grado, Italy, 2003. Available at: http://perso.ens-lyon.fr/patrick.flandrin/NSIP03 (accessed 30.08.2017).
[18] Rilling, G., Flandrin, P. One or two frequencies? The empirical mode decomposition answers // Signal Proc., IEEE Transactions. 2008. Vol. 56, No. 1. P. 85-95.
[19] Сафиуллин Н.Т. Разработка методики анализа временных рядов с помощью преобразования Хуанга —Гильберта: Дис. ... канд. техн. наук. Новосибирск, ФГОБУ ВПО, 2015. 193 с.
Safiullin, N.T. Development of the time series analysis technique using the Hilbert — Huang transform: Dis. ... kand. tekhn. nauk. Novosibirsk, FGOBU VPO, 2015. 193 p. (In Russ.)
[20] Huang, N.E., Shen, S.S.P. Hilbert — Huang transform and its applications. Singapore: World Sci. Publ. Co., 2005. 323 p.
[21] Wu, Z., Huang, N.E. Ensemble empirical mode decomposition: A noise-assisted data analysis method // Advances in Adaptive Data Analysis. 2009. Vol. 1, No. 1. P. 1-41.
[22] Torres, M.E., Colominas, M.A., Schlotthauer, G., Flandrin, P. A complete ensemble empirical mode decomposition with adaptive noise // IEEE Intern. Conf. on Acoustics, Speech and Signal Proc. (ICASSP), Prague (CZ), 2011. Prague: IEEE Publisher, 2011. P. 4144-4147.
[23] Colominas, M.A., Schlotthauer, G., Torres, M.E., Flandrin, P. Noise-assisted EMD methods in action // Advances in Adaptive Data Analysis. 2012. Vol. 4, No. 4. P. 1250025. (11 p.)
[24] Rilling, G., Flandrin, P., Goncalves, P. Empirical mode decomposition. Available at: http://perso.ens-lyon.fr/patrick.flandrin/emd.html (accessed 01.03.2012).
[25] Wu, Z., Huang, N. E. Statistical signifiance test of intrinsic mode functions // Proc. of the Royal Soc. of London A: Math., Phys. and Eng. Sci. 2004. Vol. 460, No. 2046. P. 1597-1611.
[26] Левин А.А., Таиров Э.А., Спиряев В.А. Автоколебательные пульсации давления в этаноле при захолаживании нагревателя // Теплофизика и аэромеханика. 2017. Т. 24, № 1. C. 61-72.
Levin, A.A., Tairov, E.A., Spiryaev, V.A. Self-excited pressure pulsations in ethanol under heater subcooling // Thermophys. and Aeromech. 2017. Vol. 24, No. 1. P. 61-72. (In Russ.)
Поступила в 'редакцию 28 апреля 2017 г., с доработки — 13 июня 2017 г.
72
А.А. MeBHH, B. A. CnnpaeB
Investigation of frequency-selective properties of Hilbert — Huang transform and its modifications on the example of studying the self-excited pressure oscillations
Levin, Anatoliy A., Spiryaev, Vadim A.
Melentiev Energy Systems Institute SB RAS, Irkutsk, 664033, Russia * Corresponding author: Spiryaev, Vadim A., e-mail: errolorr@gmail.com
In this paper, we consider the application of the Hilbert — Huang transform (HHT) to the analysis of non-stationary time series generated by self-excited pressure oscillations in the volume of boiling liquid. Applying the HHT allows to decompose the signal into a set of intrinsic mode functions (IMF) with its frequencies and amplitudes. However, the "mode mixing" problem appears in the study of time series of thermophysi-cal nature that are using the classical HHT. Such problem arises as the presence of oscillations of very disparate amplitude in an IMF, or the presence of very similar oscillations in different IMFs. To overcome these problems, a modified HHT method (mHHT) is proposed. It performs decomposition over an ensemble of signals plus introduction of Gaussian white noise with the subsequent averaging of the obtained IMF.
A comparison of the decomposition using HHT and mHHT by estimating the power spectral density of each IMF on the time series studied in the paper shows the absolute advantage of mHHT. In addition, since it is necessary to specify the input mHHT parameters that affect the quality of the decomposition, we made an additional study on their optimal choice. To evaluate the obtained sets of IMF corresponding to different input parameters, the concept of the qualitative set is formulated. It is based on the analysis of the frequency properties of IMF by estimating the spectral power density.
The use of mHHT with optimally selected parameters for studying the dynamics of pressure allows qualitative identification of the IMF that corresponds to the main frequencies and mechanisms of pressure oscillations. Analysis of the obtained IMF allows us to identify two stages of the self-oscillatory process: the initial stage and the stage of developed boiling.
Keywords: Hilbert — Huang transform, self-excited pressure oscillations, non-stationary time series.
Acknowledgements. The work was performed on the equipment of a unique scientific device "High-temperature circuit" ISEM SB RAS with financially support by RFBR (grant No. 15-01-01425).
Received 28 April 2017 Received in revised form 13 June 2017
© ICT SB RAS, 2017
