УДК 519
ПОСТРОЕНИЕ МОДЕЛИ СИГНАЛА КРИТИЧЕСКОЙ ЧАСТОТЫ НА ОСНОВЕ ВЕЙВЛЕТ-ПАКЕТОВ
О.В. Мандрикова, Т. С. Горева (КамчатГТУ)
Разработан подход к построению модели сигнала со сложной структурой, который базируется на разложении сигнала по ортогональным базисным функциям. В качестве базисных функций предложены семейства ортогональных вейвлет-пакетов. На основе данного подхода построена модель сигнала критической частотыfF2.
The authors describe the approach to construction of signal model with complex composition.
The basis of this approach is splitting the signal into orthogonal basic functions - waveletpackets.
Signal model of critical frequency fF2 is based on this approach.
Регистрацию данных критической частоты проводит на Камчатке ИКИР ДВО РАН. На обработке этих сигналов основываются исследования физических процессов, протекающих в ионосфере, которая является электрическим экраном Земли. Радиоволна, направленная с Земли под углом к горизонту, отразившись где-то в атмосфере, возвращается на Землю. Высота максимума электронной концентрации и сама концентрация электронов изменяются с широтой и долготой. Кроме этого, данные характеристики зависят от времени суток, сезона, солнечной активности. Для того чтобы отразить радиоволну нужной частоты, необходимо иметь в месте отражения определенную концентрацию электронов. Если она меньше необходимой величины, то радиоволна просочится через ионосферу. Критическое значение частоты отраженной радиоволны называется критической частотой. Создание адекватных моделей и автоматических методов решения поставленной задачи - серьезная проблема, трудности решения которой связаны со сложной структурой регистрируемых природных сигналов. Они являются нестационарными, включающими в себя большое разнообразие по форме и длительности локальных особенностей, связанных со многими природными явлениями и процессами (различного характера переходные процессы, аномальные эффекты, связанные с активностью Солнца, а также сопровождающие процесс подготовки сейсмических явлений и др.).
Традиционные методы, используемые для идентификации моделей временных рядов, такие как метод Бокса - Дженкинса, метод скользящего среднего, сглаживание, выделение трендов и другие, имеют, как известно, ряд ограничений, сужающих область их использования. В основе этих методов лежит предположение, что данные содержат систематическую составляющую и случайный шум (ошибку). Они включают различные способы фильтрации шума, позволяющие увидеть регулярную составляющую более отчетливо. Тем самым после применения процедуры фильтрации сигналы преобразуются в относительно гладкую кривую. В указанных методах обработки предполагается, что регулярные составляющие принадлежат к двум классам: они являются либо трендом, либо сезонной составляющей. Это влечет возникновение значительных погрешностей аппроксимации данных, накладывает ограничения на выявляемые закономерности в их поведении и ограничивает область использования.
В данной работе предложен метод построения модели сигнала критической частоты на основе конструкции вейвлет-преобразования. Эта конструкция позволяет разложить произвольный сигнал на разномасштабные составляющие и тем самым получить информацию о структуре и масштабах процесса, т. е. дает полную развертку структуры сигнала [1, 2]. Как будет показано далее, это в значительной степени повышает качество и эффективность обработки сигналов с такой сложной структурой, как сигнал критической частоты. В основе процедуры идентификации модели сигнала лежит детерминированный подход, базирующийся на разложении функции по базисам. Вейвлет-методы включают в себя широкий спектр базисных функций различной формы, в том числе ступенеобразных, пикообразных, в виде всплесков различной частоты и амплитуды. Важно отметить, что особенности в виде всплесков, ступенеобразные, пикообразные и другие, часто имеют место в сигналах регистрации различных природных параметров, в том числе и в сигнале критической частоты [2]. Кроме того, вейвлет - это единственная известная базисная функция, позволяющая реализовать аппроксимацию дискретного сигнала с наименьшей погрешностью. Разработанный авторами метод построения модели сигнала критической
частоты включает следующие этапы: разложение сигнала на составляющие на основе вейвлет-конструкции, идентификация составляющих процесса и фильтрация шумовых составляющих, оценка погрешности аппроксимации. Рассмотрим эти этапы.
I. Разложение сигнала на составляющие на основе вейвлет-конструкции
Разложение сигнала на составляющие выполнялось на основе кратномасштабного анализа (КМА) и вейвлет-пакетов (ВП) [1]. Базисные функции КМА и ВП своими сдвигами и растяжениями порождают подпространства пространства I? (К) . Пространство сигнала раскладывается на сумму подпространств, при этом сигнал единственным образом представляется в виде суммы компонент:
/ (г) = №) + Ж) +... +/ ^),
где t - время.
Компоненты сигнала / (0 имеют вид
/(О = 2с1кфЛ(0, / е 12(К),
к
где фк - базисные функции пространства I (К) .
Конструкцию КМА называют также быстрым вейвлет-преобразованием, поскольку процедура разложения заключается в реализации свертки сигнала с коэффициентами вейвлет-фильтров. Предполагая, что дискретный сигнал /[п] известен с разрешением 1, получаем, что
базисным пространством для него будет ¥0 = с1о* 2 (ф(20 t - к)), к е 2 (индекс «0» соответству-
1 (К)
ет разрешению пространства, пространство Vj имеет разрешение 2-1), функция ф называется скейлинг-функцией. Разложение сигнала /0 на две компоненты /0(^ = g_1(t) + /_^) выполняется следующим образом. Ортогональная проекция исходного сигнала /0^) на подпространство У_1 даст нам компоненту / ^), которая задается набором скалярных произведений /0(0 с функциями из ортобазиса V_1, т. е. величинами
с_-1 = ^/0(0, 2_1/2ф(2_к^ = (^т_п), 2кф(t_ 1^ = ^И,/гк+,.
Другими словами, проекция сигнала /0(1) на подпространство У_1 осуществляется путем свертки с фильтром к* (символ * обозначает транспонированный фильтр) и прореживания вдвое. В качестве компоненты g_1(t) рассматривают компоненту сигнала /0(0, ортогональную к пространству . Получают схему разложения V0 = V_l ® . Ортобазисом пространства
= \/2 2 Чк Ф^ _ к). Проекция сиг-
к
нала на подпространство Ж_1 задается набором скалярных произведений /0^) с функциями ¥ :
=^/o(t), 2_1/22 _ т^ = 2 чЛ2к+*,
что равносильно свертке с фильтром д и прореживанию вдвое.
Фильтр к* пропускает низкие частоты, тем самым осуществляя сглаживание сигнала и выделяя аппроксимирующую компоненту сигнала (эту компоненту называют также сглаженной компонентой сигнала). Фильтр д* пропускает высокие частоты, тем самым выделяя высокочастотную компоненту сигнала (эту компоненту называют также детализирующей компонентой сигнала). Описанная вычислительная процедура работает на любом масштабе и позволяет выполнить разложение сигнала на компоненты по схеме Vj = Vj_1 © . Если эту процедуру раз-
ложения выполнить для сигнала /0($) т раз, то получим следующую схему разложения:
является набор вейвлет-функций <{2 1/2 ¥| '2'_т|г, где -2-
V = ж, е ж2 е...Ф ж ФV ,
0 -1 -2 -т -т 5
где , Vj - пространства с разрешением 2 3 .
В результате исходный сигнал / будет иметь единственное представление в виде суммы
компонент:
/о (Х) = Я-1 (Х) + g-2(t) + ... + g- т (х) + /- т (Х) .
Описанная процедура разложения построена на предположении, что полезная информация о сигнале находится в низкочастотной его составляющей. В этом случае такая конструкция является эффективным инструментом выделения данной информации и представления ее в виде сглаженной компоненты сигнала. При необходимости идентификации различных типов частотно-временных структур удобным методом является конструкция вейвлет-пакетов. Куафман,
Мейер и Викерхаузер доказали, что если {0 . (23х - п)} является ортонормированным базисом
пространства ис разрешением 2-3 и к, g - пара сопряженных зеркальных фильтров, определяемых формулами 0О-1 (X) = I кп0,( 2 3 X - п ) и 03-1(?) = I gn 0 3 ( 2 ' X - п) , то семейство
п п
{00-1 (23-1 X - п), 03-1 (23-1 х - п)}п 2 является ортонормированным базисом и3 [1]. Иначе говоря, сопряженные зеркальные фильтры преобразуют ортогональный базис {0 . (23? - п)} в два ор-
тогональных семейства
{00-1 (23 1 х - п)} и {03-1 (23 1 - п)} . Пространства и- и и3-1.
по-
рожденные каждым из этих семейств нальны и и0. Ф и1 , = и . Если
ортого-положить
и3 = Ж
то, используя приведенную процедуру
имеем Ж = V. и Т • =
3 3 3
разложения этого пространства, мы получаем разбиение высокочастотной компоненты сигнала. Рекурсивное расщепление пространств на основе этой процедуры представляют в виде двоичного дерева (рис. 1), которое в вейвлет-теории называют деревом пространств вейвлет-пакетов.
Каждый узел-родитель делится на два ортогональных подпространства. На корне дерева . Применяя рекурсивное расщепление вдоль ветвей допустимого дерева, мы получаем взаимно ортогональные пространства \Ж3’1} , которые в сумме дают
: Ж0 = Ф', Ж’ .
3 3 1 1 31
Выше было сказано, что теория вейвлет-преобразования предоставляет широкий спектр базисных функций, но возможность их использования в той или иной конструкции вейвлет-разложения определяется их свойствами. Решение проблемы выбора базисной функции является нелегким и играет важную роль в задачах обработки и анализа временных последовательностей. При решении данной задачи выбор базисной вейвлет-функции осуществлялся следующим образом: на начальном этапе на основе анализа свойств вейвлетов определялось семейство базисных функций, далее выбор базиса основывался на оценке погрешности аппроксимации.
Определение семейства базисных функций. Важным свойством вейвлетов является свойство ортогональности. Оно позволяет получить компоненты сигнала, параметры которых не коррелируют между собой. Это условие повышает эффективность процедуры последующей оценки параметров модели. Помимо свойства ортогональности при выборе базиса необходимо учитывать такие важные характеристики вейвлетов, как гладкость, размер носителя, число нулевых моментов [1]. Поскольку мы решаем задачу аппроксимации, то в качестве основных критериев при выборе базиса естественным будет определить минимизацию числа аппроксимирующих слагаемых и минимизацию погрешности аппроксимации. На основе конструкции вейвлет-
разложения минимизация числа аппроксимирующих слагаемых и минимизация погрешности аппроксимации может быть достигнута выбором вейвлет-базиса, обеспечивающего как можно большее число вейвлет-коэффициентов, которые являются пренебрежимо малыми. Определяющими здесь являются такие характеристики вейвлета, как число нулевых моментов, гладкость и размер носителя. Поясним это.
1. Размер носителя. Эта характеристика вейвлета влияет на погрешность аппроксимации, особенно для конечных функций. Все известные одномерные конструкции вейвлетов приводят к базисам для 1} (К). Во многих приложениях, в частности для решения нашей задачи, интерес представляет лишь часть вещественной оси. Применяя обычные базисы вейвлетов для разложения функции, мы можем положить функцию, равной нулю вне некоторого отрезка [а, Ь], но это
порождает искусственные «скачки» на краях, находящие отражения в коэффициентах разложения. Кроме того, это неэффективно с точки зрения вычислений, т. е. чем меньше размер носителя, тем меньшую погрешность мы имеем при разложении функции. Важным здесь также является зависимость числа вейвлет-коэффициентов, не являющихся пренебрежимо малыми, от размера носителя. Если функция имеет изолированную особенность в некоторой точке Х0 и если Х0 находится внутри носителя Т к, то вейвлет-коэффициенты могут иметь большую амплитуду.
Если Т имеет носитель размера К, то при каждом масштабе 23 имеется К вейвлетов Т к, носители которых содержат Х0. Чтобы минимизировать число коэффициентов с большой амплитудой, мы должны использовать функции с наименьшим размером носителя Т .
2. Число нулевых моментов. Вейвлет Т имеет т нулевых моментов, т. е.
+да
| хк Т ({)& = 0, к = 0,т -1.
-да
Это означает, что вейвлет Т ортогонален любому многочлену степени т - 1, т. е. «пропускает» полиномы степени выше т - 1.
Размер носителя функции и число нулевых моментов априори независимы. Однако для ортогональных функций И. Добеши доказала, что если Т имеет т нулевых моментов, то его носитель имеет наименьший носитель, равный 2т - 1 [3]. Таким образом, при выборе вейвлета мы приходим к выбору между числом нулевых моментов и размером носителя. Здесь важно учитывать следующее: если сигнал имеет несколько изолированных особенностей и очень гладкий между этими особенностями, то необходимо использовать вейвлет с большим числом нулевых моментов, чтобы на малых масштабах получить большое число малых вейвлет-коэффициентов; если число особенностей нарастает, лучше уменьшить размер носителя ценой уменьшения числа нулевых моментов.
3. Гладкость вейвлета. Известным фактом является то, что число нулевых моментов и гладкость вейвлетов связаны друг с другом, но характер связи может быть различным в зависимости от вида рассматриваемого семейства базисов [3]. В работах [1, 3] показано, что, например, для сплайн-фильтров и фильтров Добеши характерно следующее свойство: гладкость вейвлета возрастает с возрастанием числа нулевых моментов. С другой стороны, Добеши в работе [3] доказала теорему, из которой следует связь нулевых моментов одной функции с регулярностью другой: если вейвлет ТеСт , то автоматически двойственный вейвлет Т имеет т нулевых моментов. Для гладких функций достижение наилучшей аппроксимации высокочастотных ее компонент обеспечивает в значительной степени большое число нулевых моментов, чем регулярность Т . Вопреки сказанному Антонини и соавторы в работе [3] на основе экспериментов показали, что в случае использования полуортогональных вейвлетов при решении задачи сжатия сигналов лучшие результаты дает более гладкий фильтр с меньшим числом нулевых моментов, чем менее гладкий фильтр с большим числом нулевых моментов. Очевидно, что полученный результат можно объяснить свойствами аппроксимируемых функций. Действительно, результат аппроксимации во многом зависит от структуры сигналов, подлежащих аппроксимации. Этот факт будет учитываться при решении данной задачи на втором этапе выбора базисного вейвлета.
Итак, мы пришли к заключению, что для выполнения выдвинутых нами требований необходимо использовать гладкие функции с наименьшим носителем и наибольшим числом нулевых моментов. Наилучшим семейством в этом смысле является семейство ортогональных вейвлетов
Добеши с компактными носителями, поскольку ортогональные вейвлеты Добеши с компактными носителями - это единственное ортогональное семейство базисных вейвлет-функций, которые имеют минимальный размер носителя при заданном числе нулевых моментов.
II. Идентификация составляющих процесса и фильтрация шумовых составляющих.
Оценка погрешности аппроксимации
Регистрируемый сигнал /[п] длиной N загрязнен добавлением шума. Шум моделируется как реализация случайного процесса Ж [п]. Тогда измеряемые данные могут быть представлены в виде
/ [п ] = /[п] + Ж [п].
Сделаем предположение, что шум Ж [п] являются белым шумом, что является вполне оправданным для большинства экспериментальных данных. Правомерность данного предположения подтверждается центральной предельной теоремой. Оценка сигнала в шуме выполняется нахождением представления, которое выделяет сигнал / из шума. Оценка вычисляется с помощью оператора, который ослабляет шум, но сохраняет сигнал. Оптимизация оператора зависит от априорной информации, имеющейся в нашем распоряжении. Например, Байесовская
процедура предполагает, что мы знаем распределение вероятности функции / , и оптимизирует оператор с целью минимизации ожидаемого риска. Такие модели часто неприменимы для сложных природных сигналов, поскольку мы не располагаем достаточной информацией для определения априорного распределения вероятности. В этом случае можно использовать минимаксную процедуру [1]. Минимаксный подход требует только знания априорного множества и гарантии того, что сигнал ему принадлежит. Для кусочно-гладких сигналов минимаксная оптимальность доказана для пороговых вейвлет-оценок [1]. Вследствие этой теории подавление шумов и восстановление сигнала было выполнено на основе пороговых оценок.
После применения конструкции вейвлет-разложения мы получаем представление сигнала в виде
Л 0) = £-1 0) + ёЛ) + ... + ё- т 0) + X т 0) ,
где ё-1 (() = 2 к^-а ({) , Л-‘ ({) = 2 с-аф-,к (') , к и ^-а - коэффициенты вейвлет-
к к
преобразования; ^_,.к - вейвлет; ф-; к - скейлинг-функция в силу ортогональности базисных функций 2 £-, (^ ) ё-к Ъ ) = Ь-;,-к и 2 ё-1 (^ )/-т (^ ) = 5-j,-т ; 5 - дельта-функция Дирака.
I I
Сигнал /0(?) аппроксимируется М векторами, адаптивно выбранными из ортонормирован-ного базиса {Т-( к, ф-( к}.=._ =_^ пространства С (Я):
ЛМ = 2(/, ^Ят .
те1м
Погрешность аппроксимации есть сумма оставшихся коэффициентов:
4М ] = ||Л - Лм||2 = 2\(/, О Г .
т£1м
Очевидно, чтобы минимизировать эту погрешность, индексы в множестве 1М должны соответствовать М векторам, имеющим скалярные произведения с наибольшей амплитудой |(/, Ят)|. Это векторы, которые наилучшим образом корректируют с / . На основе предположения, что шум Ж [п] является белым шумом, т. е. последовательностью независимых случайных величин, выделенные таким образом векторы {ят }т/ могут быть интерпретированы как наиболее характерные для / . Выделение векторов {а ) было выполнено следующим образом.
* ' те!м
Вейвлет-коэффициенты j|(f, qm}|} были упорядочены по убыванию. Обозначив
f [к] = ( f, qm^j как коэффициент ранга к : |fpr [к| > |fpr [к + 1]|, получаем аппроксимацию
M
fM = 2 f [к ]qmk .
к = 1
Данная аппроксимация сигнала далее была вычислена с применением пороговой функции:
в, (х и хесли! X >Т •
[0, если|х| < T,
где порог T таков, что f' [M + 1] < T < f' [M ]:
fM =2вт ((f • qJ )qm.
m=0
Погрешность полученной аппроксимации имеет вид
2
[M]=||/-f,,||' =2 f [к] .
к = M+1
Значение порога T естественно будет определить равным дисперсии шума W [п]. Чтобы оценить дисперсию с2 шума W [n] по данным f , мы должны подавить влияние f . Грубая оценка с2 может быть вычислена по средним значениям вейвлет-коэффициентов наименьшего масштаба [1]. Так как медиана P коэффициентов Med (а ) есть величина среднего коэф-
V Р /0<p<P
фициента ап0 ранга P, то если M — медиана абсолютного значения P независимых гауссов-
2
ских случайных переменных с нулевым средним значением и дисперсией С0 • то E {M0,6745с0 [1]. Дисперсия с2 шума W может быть оценена по медиане Mf множества f • ^_1к)} _____ пренебрежением влияния f :
Mf
с = -
0,6745
Дальнейшая степень свободы вводится адаптивным выбором базиса, зависящего от свойств сигнала.
Пусть 1'М - множество индексов М векторов рх, которые максимизируют |^/, дт)\ . Наилучшая нелинейная аппроксимация / в в х имеет вид
/м = 2( /, ят) ят.
Погрешность аппроксимации определяется как
^ ]=И|( а • чі) \ НІ АI'-£К а • «-> |г.
т&1м тєіМ
Тогда базис Р“ = {чі } 1 < т < N лучше, чем базис рт = {чЦ} 1 < т < N при аппроксимации А, если при всех М > 1
ва [М] < вт [М ].
Полученное условие лучшего базиса эквивалентно тому, что
Vм >1 £|( а , €) |’ >£|< А, €) |’.
При выборе базиса мы обеспечиваем таким способом выполнение основных критериев качества аппроксимации: минимизацию числа аппроксимирующих слагаемых и минимизацию погрешности аппроксимации. Аппроксимирующие векторы, выбранные из этого «наилучшего» базиса, наилучшим образом приспособлены для аппроксимации конкретного сигнала и описывают важные структуры сигнала, фильтруя шум. Модель сигнала в этом случае имеет вид
до = 2 2 (о+2 с-и к (о,
-,-Е/а к 1
где с_и - коэффициенты вейвлет-преобразования, соответствующие сглаженной составляющей сигнала; й_1к - детализирующие коэффициенты вейвлет-преобразования; ^ ^ - вейвлет; ф_; к - скейлинг-функция.
Decomposition Tree
’ 1 1 <з,4> (3js) <з£> (4,^0I) (4*^2) (4*13)
Рис. 2. Дерево разложения
На основе обработки данных критической частоты были получены следующие результаты:
- выделены аппроксимирующие вейвлеты Добеши 2-го порядка, наилучшим образом приспособленные для аппроксимации сигнала критической частоты, и получено дерево разложения (рис. 2);
- определена дисперсия шумовой компоненты сигнала (0,234);
- рассчитана погрешность полученной аппроксимации (78,19).
Литература
1. StephaneMallat. A Wavelet tour of signal processing / Пер. с англ. - М.: Мир, 2005.
2. Богданов В.В., Геппенер В.В., Мандрикова О.В. Моделирование временных рядов геофизических параметров на основе вейвлет-преобразования. - СПб.: Изд-во СПбГЭТУ, 2006.
3. Добеши И. Десять лекций по вейвлетам. - Ижевск: Регулярная и хаотическая динамика, 2001.