CC BY
14
ГЕОФИЗИКА
УДК 550.388.2; 520.24; 537.86
МЕТОДИКА СТАТИСТИЧЕСКОГО АНАЛИЗА ВЕЙВЛЕТ-СПЕКТРОВ ИОНОСФЕРНЫХ СИГНАЛОВ СИСТЕМЫ
GPS
В. И. Захаров, А. С. Зиенко
(.кафедра физики атмосферы) E-mail: [email protected]
Предлагается метод статистического анализа вейвлет-спектров и исследована возможность его применения для данных, полученных при радиозондировании ионосферы навигационными сигналами GPS. Показано преимущество вейвлет-анализа по сравнению со спектральным для исследуемого сигнала, в том числе и при наличии кратковременных сбоев.
Введение
В настоящее время вейвлеты стали важным математическим аппаратом во многих приложениях физики, где требуются обработка и анализ локальной структуры сложного исследуемого процесса [1-4]. В геофизике такой подход может применяться к полям различных параметров природных сред [3, 4]. В этом случае вейвлеты используются для обработки нестационарных во времени и неоднородных в пространстве сигналов, причем результат содержит не только распределение энергии в спектре сигнала, но и сведения о времени появления определенных частотных составляющих.
Базис собственных функций вейвлет-преобразо-вания обладает специальными свойствами, благодаря которым существует возможность анализировать особенности процессов, которые не могут быть выявлены с помощью традиционных преобразований Фурье и Лапласа [1-3]. Так, по сравнению с разложением Фурье вейвлеты могут представлять локальные особенности сигналов, вплоть до разрывов [1-4]. При вейвлет-преобразовании одномерных сигналов частота и время могут рассматриваться как независимые переменные, что позволяет получить сведения о временной эволюции сигнала.
1. Вей влет-спектры временных последовательностей
Прямое интегральное вейвлет-преобразование сигнала представляет собой скалярное произведение сигнала на ядро (или носитель) преобразования — функцию ф заданного типа [1-4]. Вейвлет-образ сигнала по независимым переменным а и Ь обозначим как функцию Ш(а,Ь):
W(a,b) =
s{t^ab(t) dt.
(1)
Базис функционального пространства Ь2(К) образуется путем масштабных преобразований и сдвигов вейвлета ф(1) с непрерывными значениями параметров — масштабного коэффициента а и сдвига Ь:
фаьУ) = \а\^1/2фШ-Ь)/а], а,ЬеЯ,
В качестве ядра преобразования мы использовали широко известные преобразования:
— "\¥АУЕ-вейвлет. Уравнение вейвлета с единичной нормой имеет вид
ф{1, а, Ь)
-1.786 t-b
V2
а
ехр
t-b
(2)
— МНАТ-вейвлет (Mexican hat шляпа) с уравнением
мексиканская
ф{1, а, Ь)
1.031
ехр 't-
а
t-b
ехр
t-b
(3)
— вейвлеты Морлетта до 6-го порядка включительно.
Масштабно-временная картина коэффициентов Ш(а, Ь) точно фиксирует положение структур информационного сигнала сменой знака коэффициентов Ш(а,Ь). Поэтому периодические составляющие (периодики) сигналов можно оценивать по взаиморасположению максимумов и минимумов в сечениях Ш(а, Ь) по масштабному коэффициенту (т. е. по аналогу частоты). Для самих сечений название различается у разных авторов. Так, в обзоре [3] их называют структурной или скелетной функцией (скелетоном), в [4] — скейлинг-функцией, т.е. функцией масштабов.
4? Wavelet analysis tool
□э
Start
Function wave Г MHAT С Model
Spectrum
1
0.5 0
-0,5
Wavelet Functions
I
-4 -2 0 2 4 — wave — MHAT — Morlet
Signal
I' ■' i ■' ' i'1' i1' i, 11 .11 .11,. , 1111, i' i, 11 i,
-14 -12 -10 -8 -6 -4-2 0 2 4 6 8 10 12 14 Time
Signal constructor
T1 T2 Д2/А1 T3 A3/A1
[IF [ш 10,52 [з]сГ 10,50
J]
Component Turn ON Ö РгсГ [-20" Г~ df
Noise Amplitude |l,00
-J-
Sceletons
Full Analysis
Color Scheme 1 Г RGB
<* Gray Save Graph
Параметры сигналов в смеси
Структурная функция. Сечение а = 1.40
Рис. 1. Пример анализа тестового сигнала с помощью специализированной программы. Видны заданные параметры сигналов в изучаемой смеси. В сноске показано некоторое сечение \¥(а,Ь), в окнах выводится служебная информация о найденных периодиках в данном сечении — {3.0 3.0 2.8
3.2 3.6 6.0 6.0}
Основная трудность при анализе реальных сигналов — разнообразные временные интервалы локализации информации, особенно если структура сигнала меняется во времени. Исходя из конкретной задачи, далее мы определяем оптимальный интервал преобразования, а ядро выбирается на основе решения тестовых задач, моделирующих конкретные ситуации.
2. Исследования вейвлет-спектров модельных сигналов и отработка методики статистического анализа
Для отработки методики, создания и тестирования алгоритмов автоматического выделения периодов колебаний, содержащихся в наблюдаемых данных, мы использовали специальное тестовое программное обеспечение, пример работы которого приведен на рис. 1.
В качестве тестовых процессов рассматривались различные смеси трех гармоник (с заданными периодами, амплитудами и временем включения) и аддитивного некоррелированного шума. Далее строилось вейвлет-преобразование по (1)-(3) с выбранным типом вейвлета, визуализировалось и исследовалось любое сечение а = const для функции W(a,b) преобразования.
Анализ полученного вейвлет-спектра является сложной задачей. Поэтому мы предложили и опробовали статистический подход, который заключается в исследовании структуры каждого сечения вейвлет-спектра и статистической обработке полученных результатов. Так, в сечении определяются координаты локальных максимумов и минимумов. Далее мы вводим три статистики: 1) для расстояний между минимумами, 2) для расстояний от максимума до максимума и 3) от минимума/максимума до максимума/минимума соответственно. Из полу-
Вероятность, % 35 -
30 -
25 -
20 -
15 -10 -5 -
0
Рис. 2. Распределение выделенных параметров для тестового сигнала при различных амплитудах шума. Для сравнения приведены результаты только для статистики гтпп-гтпп
чившегося ряда значений в каждой статистике для оценки периода (или полупериода — для статистики 3) выделяемой гармоники выбирается наибольшее значение. Затем производится цифровая фильтрация для удаления найденного периода. Процесс обработки повторяется до выделения минимальной периодики, и строятся различные распределения параметра «детектированный период» от частоты его появления для всех скелетонов. Под «минимальной» периодикой (и временным разрешением) здесь понимается размер ячейки гистограммы, используемой для статистической обработки. Для исследований ионосферных сигналов этот параметр естественно выбирать из физических соображений (см. п. 3), и его величина составляет не менее 3-5 мин.
На рис. 1 приведен пример анализа тестового сигнала — смеси трех компонент-синусоид: первая имеет период Т\= 6, амплитуду А\ = \ и время включения / = 0, т.е. в течение времени наблюдения присутствует чуть больше двух ее периодов. Вторая и третья гармоники действуют все время, имеют почти равные амплитуды = 0.5 и Дз = 0.52 и периоды = 3 и 7з = 1. Амплитуда аддитивного шума меняется от 0 до 3, на рисунке приведен пример, когда /4Ш„-Яе = 1. Для увеличения надежности оценки результаты обработки всех введенных статистик используются совместно. В итоге имеем распределения вероятности появления гармоник в виде, представленном на рис. 2.
Мы провели анализ влияния шума на качество определения параметров сигнала. Результаты для заданных выше параметров сигнала в функции амплитуды шума приведены на рис. 2. Получено, что предложенная статистическая методика лучше осуществляет выделение высокочастотных компонент: частота появления Р уменьшается почти линейно при увеличении периода и дополнительно падает вдвое при росте амплитуды шума в диапазоне от 0 до 3. В отсутствие шума для детектированных гармоник имеем: Я(7з) = 33%, Р(Т2) = 14%, Р{Т\) = 7%. В случае шума очень большой амплитуды /4nojSe = 2, (т.е. равной амплитуде смеси гармоник), имеем соответственно Я(7з) = 20%, Р(Т2)= 10%, Р{Т\) = 4%.
В работе использовался пороговый критерий оценки наличия структуры. Наши исследования позволяют установить такой порог появления структур в 3-4% — при автоматической обработке считается, что если выделяемая структура имеет частоту появления в W(a,b) менее пороговой, то эта периодика не рассматривается далее. Для всех уровней шума предложенная методика обеспечивает выделение заданных гармоник.
Подобные результаты предложенная методика дает при использовании wave- и mhat-вейвлетов. Сохраняются и приведенные выше зависимости вероятности появления структуры от амплитуды шума.
Амплитуда шума 7-0
2-0.5
3 - 1
4- 1.5
5-2
6- 2.5
7- 3
Применение вейвлетов Морлетта для наших целей оказалось нецелесообразным.
3. Анализ данных навигационной системы GPS
Рассмотрим использование предложенного метода вейвлет-анализа для выявления структур в ионосферных сигналах. Принято считать [5-7], что характерные частоты возбуждаемых в атмосфере и ионосфере акуетико-гравитационных волн (АГВ) лежат в диапазоне ~ 0.2 2 мГц — частоты Брен-та-Вяйсяля для типичных условий атмосферы. Использование сигналов навигационной системы GPS дает возможность выделять возмущения в ионосфере на фоне регулярных изменений на основе точных фазовых измерений [8-11], регистрируемых наблюдательной сетью IGS (International Geophysie Service). Данные имеют временное разрешение 30 с, что позволяет исследовать различные периодики сигнала, которые принято связывать с определенными процессами в ионосфере [7-8, 11-14].
Оценим необходимое время наблюдения (объем выборки). С одной стороны, оно должно превышать периоды быстрых флуктуационных процессов, с другой — необходимо, чтобы статистика исследуемого процесса была достаточной и не слишком менялась. Так получается оценка для времени наблюдения ~ 1.5^2 ч (что дает выборочную ошибку
6 -т- 8%).
Метод выделения ионосферных структур [9-14] основан на определении вариаций полного электронного содержания (ПЭС) при использовании комбинации регистрируемых фаз L\ и ¿2 в виДе
Lj = L\ -L2 = «•/ + err/, (4)
где I — наклонный по лучу зондирования ПЭС, определяемый как
1 =
1 fh
2f2
2
[LiAi — ¿2^2 + SL + const]
где « = 40.308; L]Ai и ¿2^2 приращения фазового пути на соответствующей рабочей частоте f\ и /2 с длинами волн А] и Аг [м], SL — ошибка фазовых измерений и const — неопределенность разности фаз наблюдений на разных частотах, егг/ — ошибка измерений.
Исследования [8-11] показывают, что в системе GPS ошибка в определении изменения ПЭС не превышает 1-3% (при неопределенном начальном значении этой величины). Это значение может служить оценкой величины аддитивного шума, влияние которого на эффективность работы предложенного статистического метода рассмотрено ранее.
Оценки исследуемых флуктуаций по производной процесса (4) можно построить в виде
dL[
SL[(t) = hit) - 0.5(h{t + г) + hit - г)) « -г-
dt
где величина окна цифрового фильтра г составляет от 5 мин. Мы использовали для анализа данные как для ПЭС (что обозначено индексом /), так и для регистрируемой фазы на основной частоте L\ — у нее меньше фазовый шум, а амплитуда излучаемого сигнала больше, чем у h . Полученные по (4), (5) ряды используются для спектрального и вейвлет-анализа по (1)-(3) на основе методики, представленной в п. 2.
На рис. 3 приведены результаты такого подхода. Использованы данные станции agmt за 28 сентября 2004 г. На рис. 3, а представлена двухчасовая запись анализируемого сигнала, выделенного по (4), (5) для спутника GPS N15. Рис. 3,6 показывает результаты вейвлет-преобразования изучаемого сигнала. Видно наличие разных периодик в сигнале с разной временной локализацией. На рис. 3, в приведен результат статистической обработки вей-влет-спектра по методике, изложенной в п. 2. Локальные экстремумы гистограммы, превышающие определенное ранее пороговое значение 3-4%, позволяют выделить квазигармонические составляющие в сигнале. Так выделены периодики 6 мин (соответственно частоты 2.8 мГц), 8 мин (2 мГц), 14-12 мин (1.2-1.4 мГц), 33-28 мин (0.5-0.6 мГц) и 40 мин (0.4 мГц).
На рис. 3, г приведен пример спектра мощности флуктуаций в ПЭС (сплошная линия) и в L\ (пунктир). Видны спектральные максимумы, соответствующие выделенным ранее. Трудности спектрального выделения низкочастотных компонент на данном примере: гармоника 0.4 мГц не разрешается, хотя в вейвлет-спектре она проявляется отчетливо.
Дополнительный анализ показывает, что вей-влет-преобразование возможно применять даже к сигналам, содержащим кратковременные (до 2 мин) сбои в данных. В этом случае спектральный анализ не дает результатов, а вейвлет-спектр содержит информацию о характерных периодиках, действующих в сигнале. Выбор пороговых значений при обработке в этом случае помогает игнорировать данные, связанные непосредственно со сбоем.
4. Основные результаты
Предложена методика статистической обработки вейвлет-спектров, основанная на анализе всей совокупности функций масштабов. Исследованы возможности метода при детектировании периодик на модельных сигналах. Качество выделения компонент зависит от величины шума (мешающего сигнала), но остается удовлетворительным даже при отношении сигнал/шум, равном единице. Предложен практический критерий определения порогового уровня наличия гармоники в сложном сигнале. Исследовано несколько анализирующих функций вей-влета, наилучшие результаты получены для функций типа mhat — см. (3).
Сигнал d(II3C)/dt для анализа: agmt 8:30-10:30 PRN15 2004-9-28
Карта W(a,b)
Время, ч
Выделенные вейвлет-анализом структуры d(II3C)/dt: agmt 8:30-10:30 PRN15 2004-9-28
0j j 1^20 f 40 60 .8 1.6 1.4 мГц 0.5 —0.4
80 100
Время, мин
Спектры сигналов d[II3C или Ll]/dt 2004-9-28 Время 8:30-10:30 PRN15 Фильтр 5.0 мин
Рис. 3. Пример сравнения методик обработки сигналов системы GPS: а — экспериментальный сигнал, 28.09.2004, 08:30-10:30 UTC, GPS PRN15, станция AGMT; б - карта коэффициентов W(a,b); в — результат статистической обработки вейвлет-спектра; г — результаты спектрального анализа
исследуемых сигналов (ПЭС — сплошная кривая)
Методика опробована на реальных ионосферных сигналах, полученных при обработке данных GPS. Сравнение спектральной и вейвлет-обработки выявило преимущества последней в низкочастотной части спектра. Анализ показывает, что вейвлет-об-работка эффективна даже для сигналов, содержащих кратковременные сбои.
Изложенные в статье методики реализованы в разработанном авторами комплексе CRASS (Complex Region Analysis Satellite Signals) GPS — ком-
плексного регионального анализа спутниковых сигналов GPS.
Работа выполнена при частичной финансовой поддержке РФФИ (гранты 06-05-64988а и 04-0564671а).
Литература
1. Добеши И. Десять лекций по вейвлетам. М., 2001.
2. Чуй К. Введение в вейвлеты. М., 2001.
3. Астафьева Н.М. // УФН. 1996. 166, № 11. С. 1145.
4. Дремин И.М., Иванов О.В., Нечитайло В.А. 11 УФН. 2001. 171, № 5. С. 465.
5. Госсард Э., Хук У. Волны в атмосфере. М., 1978.
6. HockeK., Schlegel К. // Ann. Geophys. 1996. 14, N 5. P. 917.
7. Ахмедов P.P., Куницын B.E. 11 Вестн. Моск. ун-та. Физ. Астрон. 2003. № 3. С. 38 (Moscow University Phys. Bull. 2003. N 3. P. 49).
8. Hoffmann- Wellenhof В., Lichtenegger H., Collins J. Global Positioning System: Theory and Practice. N.-Y., 1992.
9. Davies K., Hartmann G.K. 11 Radio Sci. 1997. 32. P. 1695.
10. Fitzgerald T.J. 11 J. Atm. and Solar-Terr. Phys. 1997. 59. P. 829.
11. Куницын В.E., Андреева Е.С., Кожарин М.А. и др. // Вестн. Моск. ун-та. Физ. Астрон. 2005. № 1. С. 74 (Moscow University Phys. Bull. 2005. N 1. P. 94).
12. Pi X., Mannucci A.J., Lindgwister U.J., Ho CM. 11 Geophys. Res. Lett. 1997. 24. P. 2283.
13. Afraimovich E.L., Palamartchouk K.S., Perevalo-va N.P. 11 J. Atm. Terr. Phys. 1998. 60, N 12. P. 1205.
14. Ахмедов P.P., Куницын B.E. 11 Геомагнетизм и аэрономия. 2004. 44", № l.C. 1.
Поступила в редакцию 20.03.06
