CC BY
0Всероссийская открытая научная конференция «Современные проблемы дистанционного зондирования, радиолокации, распространения и дифракции волн» - Муром 2024
УДК 51-73 DOI: 10.24412/2304-0297-2024-1-478-487
Применение метода эмпирической модовой декомпозиции для обработки сигналов сложной формы
И.А. Насыров
Казанский федеральный университет, Институт физики
420008, г. Казань, ул. Кремлёвская, 16а.
E-mail: Igor.Nasyrov@kpfu.ru
В статье описывается метод эмпирического модового разложения, позволяющий обрабатывать и анализировать нелинейные и нестационарные сигналы. Рассмотрено применение данного метода для удаления общих трендов и фильтрации временных последовательностей экспериментальных данных, полученных при помощи ГНСС-приёмников и фотометров во время проведения нагревных экспериментов на стенде СУРА. Ключевые слова: нелинейные нестационарные сигналы, метод эмпирической модовой декомпозиции, полное электронной содержание, ПЭС, фотометрическая кривая
Application of the empirical mode decomposition method for processing complex shape signals
I.A. Nasyrov
Institute of Physics, Kazan Federal University.
The method of empirical mode decomposition that allows processing and analysing nonlinear and nonstationary signals are described in the paper. The application of this method for removal of general trends and filtering of time sequences of experimental data obtained by GNSS receivers and photometers during heating experiments at the SURA facility is considered. Keywords: nonlinear and nonstationary signals, method of empirical mode decomposition, total electron content, TEC, photometrical curve.
Введение
Большинство естественных материальных процессов, реальных физических систем и соответствующих этим процессам и системам данных в той или иной мере являются нелинейными и нестационарными, и при анализе данных используются определенные упрощения, особенно в отношении априорно устанавливаемого базиса преобразования данных в новые, удобные для обработки и анализа метрические пространства [1]. И только в последние десятилетия начали активно развиваться методы анализа нелинейных, но стационарных и детерминированных систем, и линейных, но нестационарных данных (вейвлет анализ [2], распределение Вигнера-Виля (см. например [3]) и др.).
Таким образом, в природе мы чаще всего имеем дело с нелинейными и нестационарными сигналами, для анализа которых необходим адаптивный базис, получаемый в ходе определенного метода из самого сигнала. Такой метод был предложен Норденом Хуангом в 1995 году и назван методом эмпирической модовой декомпозиции (ЭМД). В 1998 году метод был расширен за счет преобразования Гильберта и обобщен для анализа любых временных данных [4-6]. Этот метод получил название - преобразование Гильберта-Хуанга (Hilbert-Huang Transform, HHT), он часто используется в исследованиях изменения климата, океанических волн, при анализе спутниковых, геофизических, метеорологических и биомедицинских данных и т. д. [7-
9]. Этот метод является чисто эмпирическим и не нуждается в априорной информации, что делает его высокоадаптивным к различным задачам [10]. В ряде работ показано, что метод эмпирической модовой декомпозиции превосходит вейвлет-анализ по частотно-временному разрешению [11].
Описание метода эмпирической модовой декомпозиции сигналов
Эмпирическая модовая декомпозиция (Empirical Mode Decomposition, EMD) - метод разложения сигналов на функции, которые получили название внутренних или «эмпирических мод». Метод представляет собой адаптивную итерационную вычислительную процедуру разложения исходных данных (непрерывных или дискретных сигналов) на эмпирические моды или внутренние колебания.
У каждого сигнала имеются локальные экстремумы: чередующиеся локальные максимумы и локальные минимумы с произвольным расположением по координатам (независимым переменным) сигналов. По этим экстремумам с использованием методов аппроксимации можно построить две огибающие сигналов: нижнюю - построенную по точкам локальных минимумов, и верхнюю - построенную по точкам локальных максимумов, а также функцию «среднего значения огибающих», которой отвечает срединная линия, расположенная в точности между нижней и верхней огибающими.
Модовая декомпозиция сигналов основана на предположении, что любые данные состоят из различных внутренних колебаний (intrinsic mode functions, IMF). В любой момент времени данные могут иметь множество сосуществующих внутренних колебаний - собственных модовых функций (СМФ). Каждое колебание, линейное или нелинейное, представляет собой модовую функцию, которая имеет экстремумы и нулевые пересечения. Кроме того, колебания в определенной степени «симметричны» относительно локального среднего значения. Конечные сложные данные образуются суммой СМФ, наложенных на региональный тренд сигнала.
Необходимо, чтобы каждая СМФ соответствовала следующим утверждениям:
1) разность между общим количеством экстремальных точек и количеством пересечений данной функции с осью абсцисс не должна превышать единицы;
2) в любой точке функции среднее значение между их верхней и нижней огибающими, полученными путем аппроксимации локальных максимумов и минимумов соответственно, не должно быть отличным от нуля.
СМФ представляет собой колебательный режим, но вместо постоянной амплитуды и частоты, как в простой гармонике, у СМФ могут быть переменная амплитуда и частота, как функции независимой переменной (времени, координаты, и пр.). Первое свойство гарантирует, что локальные максимумы функции всегда положительны, локальные минимумы соответственно отрицательны, а между ними всегда имеют место пересечения нулевой линии. Второе свойство гарантирует, что мгновенные частоты функции не будут иметь нежелательных флуктуаций, являющихся результатом асимметричной формы волны. Любой произвольный сигнал, изначально содержащие произвольную последовательность локальных экстремумов (минимум 2), можно разделить на семейство СМФ и остаточный тренд. Если данные лишены экстремумов, но содержат точки перегиба («скрытые» экстремумы наложения модовых функций и крутых трендов), то для открытия экстремумов может использоваться дифференцирование сигнала.
Допустим, что имеется произвольный сигнал y(t). Сущность метода ЭМД заключается в последовательном вычислении функций эмпирических мод Cj(t) и остатков rj(t) = Tj-1(t) - Cj(t), где j = 1, 2, 3, ..., n при ro = y(t). Результатом разложения будет представление сигнала в виде суммы модовых функций и конечного остатка:
у(1) = + Гп(1),
¡=1
где п - количество эмпирических мод, которое устанавливается в ходе вычислений.
Алгоритм эмпирической декомпозиции сигнала складывается из следующих операций его преобразования:
1) находим в сигнале у(к) положение всех локальных экстремумов, максимумов и минимумов процесса (номера точек к^.ехг экстремумов), и значения у^ких) в этих точках. Между этими экстремумами сосредоточена вся информация сигнала. Группируем раздельно для максимумов и для минимумов массивы координат кег и соответствующих им амплитудных значений у(к1.ехл). Число строк в массивах максимумов и минимумов не должно отличаться более чем на 1;
2) кубическим сплайном (или каким-либо другим методом) вычисляем верхнюю иг(к) и нижнюю иь(к) огибающие процесса соответственно, по максимумам и минимумам. Определяем функцию средних значений ш\(к) между огибающими:
... щ(к) + иь(к) т1(к) =---.
Разность между сигналом у(к) и функцией ш\(к) дает нам первую компоненту отсеивания - функцию М(к), которая является первым приближением к первой функции СМФ:
= у(к) -т1(к);
3) повторяем операции 1 и 2, принимая вместо у(к) функцию М(к), и находим второе приближение к первой функции СМФ - функцию к2(к):
Ь.2(к) = к1(к) -т2(к);
последующие итерации выполняются аналогично. Алгоритм итераций нахождения первой функции СМФ:
к1(к) = к1-1(к) -т1(к).
По мере увеличения количества итераций функция даг(к) стремится к нулевому значению, а функция Иг(к) - к неизменяемой форме. С учетом этого, естественным критерием останова итераций является задание определенного предела (с>) по нормализованной квадратичной разности между двумя последовательными операциями приближения.
Как правило, для выполнения качественного отсеивания модовых функций достаточно 6-8 итераций. Слишком строгий критерий останова (с>) может завышать количество СМФ и создавать компоненты, не несущие какой-либо полезной информации. С другой стороны, при слабом критерии возможно отсеивание СМФ, не полностью удовлетворяющих свойствам модовых функций. Останов итераций по нормализованной квадратичной разности исторически был первым по применению. В 2003 г. Quek [12] предложил другой, более эффективный критерий переменной меры, определенной как
1к1к1-1(к)-к1(к)12
8 =
Последнее значение Иг(к) итераций принимается за наиболее высокочастотную функцию С1(к) = Иг(к) семейства СМФ, которая непосредственно входит в состав
исходного сигнала у(Ь). Это позволяет вычесть из состава сигнала и оставить в нем более низкочастотные составляющие:
ч(к) = у(к) - ci(k).
Функция п(£) обрабатывается как новые данные по аналогичной методике с нахождением второй функции СМФ - C2(k), после чего процесс продолжается:
Таким образом, достигается декомпозиция сигнала в п - эмпирическом приближении:
Остановка декомпозиции сигнала должна происходить при максимальном «выпрямлении» остатка, т.е. превращения его в тренд сигнала по интервалу задания с числом экстремумов не более 2-3. Даже для данных с нулевым средним значением конечный остаток может отличаться от нуля. Чтобы применять метод ЭМД, центрирования данных не требуется, метод нуждается только в локализациях экстремумов. Нулевая линия для каждого компонента декомпозиции формируется процессом отсеивания. Извлеченные СМФ локально симметричны и имеют физически значимые функции мгновенных частот. Каждая СМФ содержит более низкие частотные составляющие, чем извлеченная перед ней.
Далее рассмотрим два примера применения метода эмпирической модовой декомпозиции для удаления общего тренда в сигналах сложной формы.
Удаление общего тренда из данных измерений ПЭС
Измерение полного электронного содержания (ПЭС) по фазовым данным приема сигналов глобальных навигационных спутниковых систем (ГНСС) в настоящее время является одним из традиционных методов решения различных геофизических задач [13]. Измерения ПЭС используются при изучении источников возмущений ионосферной плазмы, таких как солнечный терминатор, солнечные вспышки, геомагнитные бури, землетрясения, падение метеоров, запуски ракет, воздействие мощного наземного радиоизлучения и другие [14-17]. Изучение различных ионосферных неоднородностей является сегодня актуальной задачей, поскольку эти неоднородности могут влиять на распространение радиоволн и существенно влиять на надежность и помехоустойчивость наземных и космических радиосистем. Однако для получения информации о величине и масштабе регистрируемых неоднородностей, наряду с проведением измерений абсолютных значений ПЭС, необходимо выделять вариации ПЭС, вызванные различными ионосферными процессами. Для этого, прежде всего, необходимо удалить из данных ПЭС тренд, который зависит от изменения угла возвышения спутника и, соответственно, связан с расстоянием от этого спутника до ГНСС приемника, а затем выполнить процедуру цифровой фильтрации. Для удаления тренда чаще всего используются такие методы, как скользящее среднее и вычитание аппроксимирующих полиномов. Для выделения вариаций ПЭС используются такие методы цифровой фильтрации, как скользящее среднее, фильтры Баттерворта различных порядков и другие.
Продемонстрируем применение метода эмпирической модовой декомпозиции для удаления общего тренда в записи ПЭС (рис. 1 б), полученной на стенде СУРА при пролёте навигационного спутника GPS G08 через главный лепесток диаграммы
Г2(к) = 1 (Ю -С2(к), и т.д.
направленности антенной системы (ДН) во время нагревных экспериментов 15 марта 2010 г. (рис. 1 а). По оси ординат на рис. 1 б) вариации ПЭС указаны в общепринятой системе единиц ТБСИ (1 ТБСИ = 1016 эл/м2).
ПЭС, TECU
Рис. 1. Пролёт навигационного спутника GPS G08 через ДН стенда СУРА 15.03.2010 во
время проведения нагревнго эксперимента: а) траектория движения проекции подспуниковой точки. Маркерами обозначены двухминутные временные интервалы, овал - проекция ДН стенда СУРА на высоте движения спутника; б) наклонный ПЭС, полученный по фазовым измерениям на ГНСС станции, расположенной на стенде СУРА.
Вертикальными штриховыми линиями показан временной интервал пролёта спутника
через ДН стенда СУРА.
Для реализации метода эмпирической модовой декомпозиции была разработан специальный программный пакт. Применение данного программного пакета к записям ПЭС, полученным при помощи ГНСС приёмника, позволило эффективно удалить остаточный тренд из наклонных данных ПЭС (красная кривая на рис. 2 а).
В дополнение к удалению остаточного тренда, использование данного метода позволило выделить из данных ПЭС флуктуации различных масштабов, содержащиеся в сигнале в виде отдельных СМФ, упорядоченных по частоте и различных по своей магнитуде. В данном конкретном случае пролёта навигационного спутника GPS G08 через ДН стенда СУРА во время проведения нагревных экспериментов было выделено десять отдельных мод (рис. 2 б), каждая из которых, по всей видимости, соответствует отдельному физическому процессу, оказывающему влияние на ионосферу. Ввиду интегрального характера измерений ПЭС, изучение отдельных СМФ, а также физических процессов, ответственных за их формирование и вклад в общую картину ПЭС, требует отдельного детального исследования и выходит за рамки настоящей работы.
Как упоминалось ранее, получение остаточного тренда с помощью метода эмпирической модовой декомпозиции является результатом выделения СМФ из данных ПЭС, которые в сумме представляют вариации ПЭС после удаления тренда (рис. 2 в).
Рис. 2. Применение метода ЭМД к цифровой фильтрации данных наклонного ПЭС: а) - зависимость наклона ПЭС от времени показана чёрной линией, а остаточный тренд, полученный методом ЭМД, - красной линией; б) - все СМФ, отобранные с помощью метода ЭМД; в) - вариации ПЭС после удаления тренда.
Фильтрация фотометрических данных, полученных на стенде СУРА
Одним из эффектов воздействия мощных высокочастотных электромагнитных полей на ионосферную плазму является увеличение интенсивности оптического свечения ионосферы. Под влиянием радиоволны накачки с обыкновенной поляризацией в области её отражения генерируются плазменные волны, которые могут эффективно ускорять электроны. Если ускоренные электроны обладают достаточной энергией, то, сталкиваясь с атомами ионосферного газа, они возбуждают определённые энергетические уровни этих атомов. Оптическое свечение возникает при переходе электронов с возбуждённых уровней на более низкие энергетические уровни. Основные измерения оптического свечения проводят в красной (Л = 630 нм) и зелёной (Л = 557,7 нм) линиях атомарного кислорода, связанных с переходами электронов с уровня О(^) в основное состояние 0(3Р) и с уровня О(^) на уровень 0(1Б) соответственно. Регистрация искусственного оптического свечения ионосферы может применяться для решения целого ряда задач: для оценки функции распределения и концентрации ускоренных электронов, исследования дрейфовых движений в возмущённой области, изучения характеристик пятна свечения в зависимости от различных ионосферных и аппаратных условий (ориентации диаграммы направленности, характеристик образуемой возбуждённой области и т. п.), регистрации искусственной ионизации в возбуждённой области, исследования крупномасштабной структуры возмущённой области ионосферы [18 - 23].
Регистрация оптического свечения ионосферы осуществляется, в основном, при помощи ПЗС камер, которые снимают портреты ночного неба, а также при помощи высокочувствительных электрофотометров, которые записывают величину светового потока. В обоих случаях регистрирующие приборы снабжаются узкополосными, как правило интерференционными, светофильтрами, настроенными на определённую линию оптического диапазона.
На рис. 3 а представлена фотометрическая кривая в зелёной линии оптического спектра (Л = 55,7 нм), полученная при помощи высокочувствительного фотометра, установленного в непосредственной близости от стенда СУРА. Измерения проводились 10.09.2023, воздействие на ионосферу осуществлялось мощной радиоволной О-поляризации, частота волны накачки варьировалась согласно программе проведения измерений от 5,26 до 5,38 МГц в режиме 2,5 мин - нагрев, 3,5 мин - пауза.
На поведение фотометрической кривой, представленной на рис. 3 а, влияют одновременно несколько процессов. Во-первых, наблюдается экспоненциальное понижение общего фона свечения неба, связанное с заходом Солнца. Во-вторых, в фотометрической записи присутствуют вариации свечения в зелёной линии, связанные с изменениями фона естественной эмиссии верхней атмосферы на высотах ~100 км. В-третьих, в записи присутствует интегрированный световой поток от звёзд, попавших в поле зрение фотометра во время проведения измерений. Кроме этого, надо сказать, что тепловой шум фотоэлектронного умножителя (ФЭУ), на основе которого построен фотометр, имеет пуассоновское распределение и его величина зависит от величины регистрируемого светового потока. Таким образом, искусственное свечение ионосферы, стимулированное мощным радиоизлучением стенда СУРА накладывается на все вышеперечисленные естественные вариации. Здесь необходимо заметить, что все естественные вариации вызывают увеличение светового потока на сотни и даже тысячи Релей, а увеличение светового потока, связанное со стимулированным свечением ионосферы в зеленой линии оптического спектра (Л = 557,7 нм) может составлять лишь несколько единиц Релей. Это означает, что при регистрации стимулированного свечения ионосферы реализуется низкое соотношение сигнал/шум, так как полезным сигналом в данном случае является уровень стимулированного свечения.
На рис. 3 б представлено разложение фотометрической кривой (рис. 3 а) на СМФ.
Рис 3. Фотометрическая кривая свечения атмосферы в злёной линии атомарного
кислорода (Л = 55,7 нм), записанная 10.09.2023 во время проведения нагревных экспериментов на стенде СУРА (а). Разложение фотометрической кривой на СМФ при помощи метода эмпирической модовой декомпозиции (б)
Как можно видеть на рис. 3 б фотометрическая кривая при помощи метода ЭМД разложилась на одиннадцать СМФ плюс остаток. Очевидно, что СМФ №1 и СМФ №2 характеризуют тепловой пуассонов шум самого ФЭУ, поэтому их можно исключить из дальнейшего рассмотрения. Остаток, а также СМФ №№ 8-11 описывают достаточно сложный общий тренд фотометрической кривой.
Таким образом, для дальнейшей обработки фотометрических данных необходимо взять моды с третьей по восьмую. Просуммируем эти кривые, а затем при помощи
нахождения локальных отрицательных минимумов вернём все фотометрические отсчёты в положительную полуплоскость. Результаты такой обработки сигнала представлены на рис. 4.
А, отн.ед.
1 1 1 1 1 1 1 i ill II u H f 1 1 1 1 , IJ 1 1 1
1 Jf |i n 1 ЩЩЩ HW "ШН IK ll f, ч kW W W W If * W Щ 1 H
16:50 17:02 17:14 17:26 17:38 17:50 18:02 18:14 18:26 18:38 18:50 19:02 19:14 19:26 19:38 19:50 20:02 20:14 UTC
Рис. 4. Зелёная линия - результаты сложения СМФ №№ 3-8, представленных на рис. 3 б; чёрная линия - моменты включения выключения мощного радиоизлучения стенда СУРА
Как можно видеть из рис. 4 периоды увеличения искусственного свечения ионосферы в зелёной линии оптического спектра (Л = 55,7 нм), стимулированного мощным радиоизлучением стенда СУРА очень хорошо выделяются. при помощи метода эмпирической модовой декомпозиции сигналов сложной формы.
Заключение
Таким образом, в заключение можно сказать, что метод эмпирического модового разложения сигналов сложной формы может с успехом применяться при цифровой обработке как фотометрических данных, так и данных по вариациям ПЭС, полученным при помощи ГНСС приёмников.
Работа выполнена при поддержке гранта РНФ №23-27-00323.
Литература
1. Марпл-мл С.Л. Цифровой спектральный анализ и его приложения: Пер. с англ. - М.: Мир, 1990. - 584 с.
2. Малла С. Вэйвлеты в обработке сигналов: Пер. с англ. - М.: Мир, 2005. - 671 с.
3. Шкелёв Е.И. Частотно-временное разрешение автономного и взаимного преобразований Вигнера-Виля // Известия вузов. Радиофизика. - 2015. - Т. 58, № 1. - С. 50-61.
4. Huang N. E., Shen Z., Long S. R., Wu M. C., Shih H. H., Zheng Q., Yen N. C., Tung С. C., and Liu H. H. The empirical mode decomposition and the Hilbert spectrum for nonlinear and nonstationary time series analysis // Proc. R. Soc. Lond. A, 454, 1971, 1998, pp. 903-995, doi: 10.1098/rspa.1998.0193.
5. Huang N. E., Shen S.S.P. The Hilbert-Huang transform and its applications: World Scientific Publishing Co. Pte. Ltd., 2005, 311 p, doi: 10.1142/5862.
6. Huang N. E., Wu Z. A review on Hilbert-Huang transform: Method and its applications to geophysical studies // Rev. Geophys., 46(2), 2008, p. RG2006. doi: 10.1029/2007RG000228.
7. Coughlin K. T., Tung K. K. 11-year solar cycle in the stratosphere extracted by the empirical mode decomposition method // Adv. Space Res., 34(2), 2004, pp. 323-329, doi: 10.1016/j.asr.2003.02.045.
8. Souza Neto E. P., Custaud M. A., Cejka J. C. et al. Assessment of cardiovascular autonomic control by the empirical mode decomposition //Method. Inform. Med., 43, 1, 2004, pp. 60-65, doi: 10.1055/s-0038-1633836.
9. Wu Z., Huang N. E. A study of the characteristics of white noise using the empirical mode decomposition method // Proc. R. Soc. Lond. Ser. A., 460, 2046, 2004, pp. 1597-1611, doi: 10.1098/rspa.2003.1221.
10. Flandrin P., Rilling G., Goncalves P. Empirical mode decompositions as a filter bank // IEEE Signal Process. Lett., 11, 2, 2004, P. 112-114, doi: 10.1109/LSP.2003.821662.
11. Huang N. E., Shen Z., Long S. R. A new view of nonlinear water waves: the Hilbert spectrum // Annu. Rev. Fluid Mech., 31, 1, 1999, pp. 417-457, doi: 10.1146/annurev.fluid.31.1.417.
12. Quek, S., Tua, P., and Wang, Q. Detecting anomalies in beams and plate based on the Hilbert-Huang transform of real signals // Smart Materials and Structures 12, 2003, pp. 447460.
13. Bust G. S., Mitchell C. N. History, current state, and future directions of ionospheric imaging // Rev. Geophys., 46, 1, 2008, pp. RG1003, doi: 10.1029/2006RG000212
14. Afraimovich E. L., Edemskiy I. K., Voeykov S. V., Yasyukevich Y. V., Zhivetiev I. V. The first gps-tec imaging of the space structure of ms wave packets excited by the solar terminator // Ann. Geophys., 27, 4, 2009, p. 1521-1525, doi: 10.5194/angeo-27-1521-2009
15. Tsugawa T., Otsuka Y., Coster A. J., Saito A. Medium-scale traveling ionospheric disturbances detected with dense and wide tec maps over north America // Geophysical Research Letters, 34, 22, 2007, p. L22101, doi: 10.1029/2007GL031663.
16. Tsugawa T., Nishioka M., Ishii M. et al. Total electron content observations by dense regional and worldwide international networks of gnss // Journal of Disaster Research, 13, 3, 2018, pp. 535-545, doi: 10.20965/jdr.2018.p0535.
17. Kunitsyn V. E., Padokhin A. M., Vasiliev A. E. et al. Study of GNSS-measured ionospheric total electron content variations generated by powerful HF-heating // Advances in Space Research, 47, 10, 2011, pp. 1743-1749, doi: 10.1016/j.asr.2010.03.031.
18. Grach S. M., Nasyrov I. A., Kogogin D. A. et al. Mutual allocation of the artificial airglow patches and large-scale irregularities in the hf-pumped ionosphere // Geophysical Research Letters, 45, 12, 2018, pp. 12,749-12,756, doi: 10.1029/2018GL080571
19. Gustavsson B., Sergienko T., Kosch M. J., et al. // Ann. Geophys. 2005. V. 23. P. 1 747.
20. Grach S.M., Sergeev E.N., Nasyrov A.M., et al. // Adv. Space Res. 2004. V. 34, No. 11. P. 2 422.
21. Grach S.M., Kosch M. J., Yashnov V.A., et al. // Ann. Geophys. 2007. V. 25. P. 689.
22. Kosch M. J., Pedersen T., Rietveld M.T., et al. // Adv. Space Res. 2007. V. 40. No. 3. P. 365.
23. Грач С.М., Клименко В. В., Шиндин А. В. и др. // Изв. вузов. Радиофизика. 2012. Т. 55, № 1-2. С. 36.
