Применение вейвлет-анализа к задачам исследования загрязнения окружающей среды Текст научной статьи по специальности «Математика»

УДК 517. 444 : 551 . 510. 42

В.Г.Захаров, П.Г.Фрик

ПРИМЕНЕНИЕ ВЕЙВЛЕТ-Д11АЖЗА К ЗАДАЧАМ ИССЛЕДОВАНИЯ ЗАГРЯЗНЕНИЙ ОКРУЖАЮЩЕЙ СРЕДЫ

Abstract

Recently a new type of functional basis has been introduced in analysis, wavelets, which permit an analysis both in physical space and in scale space. A short revieu) of main properties of пхме let-transform is done. The method is applied to data of atmosphere pollution. The choice of best basis to obtain a minimal number of basis functions is discussed.

Среда задач по метеорологическим аспектам загрязнения атмосферы наибольший интерес представляют вопросы краткосрочного и долгосрочного прогноза состояния окружающей среды. Задачи прогноза требуют построения модели сложной нелинейной системы, зависящей как от многочисленных метеорологических, так и техногенных факторов, В этой связи большое значение приобретают исследования закономерностей пространственно-временного распределения концентраций примесей, что стимулирует в свою очередь использование разнообразных методов для анализа аэрометрических данных. Среди прочих используется и спектральный анализ, но его возможности в таких задачах весьма скромны и сводятся к выявлению достаточно очевидных периодических составляющих, соответствующих суточным, недельным ш годовым циклам сгз. На практике встречаются ситуации, когда чередуются различные типы суточной изменчивости концентраций и необходимо разделить ооответстующие временные интервалы, либо требуется выделить разовые выбросы, происходящие на фоне периодических и квазипериодических сигналов. При построении математической модели процесса заманчиво иметь базис, подстраивающийся под текущую структуру сигнала.

Альтернативой фурье-анализу становится возникший совсем

недавно математический аппарат, названный вейвлет-анализом и уже успешно конкурирующий с анализом Фурье в таких' областях, как обработка и синтез сигналов и изображений, распознавание образов, исследование функциональных операторов и сложных многомерных полей, изучение турбулентности и др. В основе вейвлет-анализа лежит представление исследуемой функции по базису, каждая функция которого характеризует как определенную временную с пространственную:» частоту, так и место ее локализации во времени с пространствен. Это достигается за с чет выбор;) специальных функций - вейвлетов С wavelet. - небольшая волна, рябьз, локализованных как в физическом пространстве, -:к и в пространстве частот и подучаемых друг из друга путем масштабного преобразования и сдвига.

Идеи использования функций, попадающих под определение г>ейвлета, высказывались в разное время разными авторами и восходят по сути к работам Хаара 1900 года, но целенаправленное развитие теории вейвлетов началось после работ А.Гроссмана и Ж. Морле г вj с 1984г.3, стимулированных проблемами обработки результатов сейсмических наблюдений. К настоящему времени число статей по вейвлетам превысило 1000, имеются обзоры с 7з и монографии tблгз. К сожалей. вейвлет-анализ мало известен в нашей стране, хотя после работы Зимина 1881г. езз, подложившего для описания турбулентности иерархический базис сфактически веивлет-базисо, еще до появления теории вейвлетов была выполнена большая серия работ, в которых система разномасштабных функций использовалась не только для анализа, но и для моделирования нелинейных процессов в развитой турбулентности, с Обзор этих работ CM. В [4.813

Цель данной статьи состоит в том, чтобы, кратко изложив основные идеи вейвлет-анализа, рассмотреть возможности его применения к задачам исследования аэрометрических данных. Особое внимание уделяется проблеме выбора базиса, оптимального с точки зрения описания системы с помощью минимального числа функций.

ВЕйЬЛЕТЫ. .Мая вейвлет-функция данного семейства

подучается из единственной "материнской" функции у путем сжатия

< растяжения'1 и сдвига:

-1 -^2 -1

V/ = ^ уЛ х “ Ь Л > <. 1 >

га,Ь ¥

где пор::'"-:/гр а « к* - масштабный множитель, отвечающий за

ширину ве>,л^:т,а и ь .= к - параметр сдвига, соответствующий

положению вейвлета. Таким образом, все функции базиса должны быть самоподобяы и. частности, имеют постоянное число

осцилящй. Благодаря это**.}- условию, вейвлет-преобразование дает хорошее пространственное разрешение при малых масштабах и хорошее частотное разрешение при больших.

С Формальной точки зрения семейство вейвлетов можеч

строиться на основе любой функции у. удовлетворяющей условию

сс

11^С х5е2.< — О. С2>

-00

На практике вейвлет должен также быть хорошо локализован как в физическом пространстве, так и в пространстве частот. Дя» этого достаточно, чтобы функция у- была сконцентрированна в некоторой конечной области и была достаточно гладкой. Специфика локализации различных функций в физическом и Фурье-пр(«-.гранстьах показана на рис. 1.

Рис. 1 са:> Вейвлет мексиканская шляпа; сб:> функция лиглвуда--Пели. функция Хаара; с г> вейвлет Мейера.

Необходимо отметить особенности вейвлет-преобразования, выгодно отличающие его от преобразования Фурье. Последнее распределяет сингулярные возмущения по всему спектру частот, тем самым делая очень трудным или даже невозможным изучение .локальных свойств сигнала из его фурье-образа. Если эти возмущения имеют случайный характер. то их практически невозможно отфильтровать. Вейвлет-преобразование сохраняет локальность подставления сигнала и позволяет локально его и восстановить. Возможно реконструировать только часть сигнала или выделить вклад определенного масштаба. Если вейвлет-коэффициенты подвержены случайным ошибкам. они будут действовать на реконструируемый сигнал локально вблизи положения возмущения. Свойства вейвлет-анализа обусловливают’ его эффективность для выделения в сигнале непериодических и квазипериодических

< труктур различно!о временного и/или пространственного масштаба.

КГЯРЕРЫВНОЕ И ДИСКРЕТНОЕ ПРЕОБРАЗОВАНИЕ. Непрерывное вейвлет-преобразование wC а, Ю фуНКЦИИ fCx}. х е R еСТЬ

jd

Г -1/3 * г'.. - ь!

v*. л. bJ - С V 1 J - , л -----IlCxDdx, С 33

j . а J

—on

Условие сг:> гарантирует обратимость преобразования: должна существовать по крайней мере одна формула для точной реконструкции сигнала из его вейвлет-образа. В качестве примера функции, используемой для непрерывного преобразования, приведем вейвлет Марра с мексиканская шляпа:», который является второй производной функции Гаусса: y-cx-j = о -х2зехР{-х2/г> с рис. о.

Характерной особенностью непрерывного преобразования является его избыточность, вытекающая из того, что пространство функций одного переменного отображается в .двухмерное пространство • масштаб-координата'-», и прояавдющаяся в наличии корреляции между веивлет-коэффициентами. Непрерывное преобразование эффективно применяется для анализа многомерных полей, обработки изображений и звуковых сигналов. Избыточность делает его устойчивым к возможным случайным возмущениям, однако корреляции

между коэффициентами могут привести к неверной интерпретации результатов преобразования исследуемого сигнала. Кроме того, его реализация требует значительных вычислительных ресурсов.

Ь случае дискретного вейвлет-преобразования масштабный параметр а и параметр сдвига ь принимают дискретные значения:

а = лт, Ь = пЬ ат, С*,п 6 л; а >1 ; Ь >0 } . ТОГДЭ о о о о о ^

-тъ/2 -го , . гп~ - -туг ^ -~тг» .

ш С х) = а и* а с. х - пЬ а 3.) = а ц<. а х - пЬ .

™,п о о оа о о о

Дискретное вейвлет-преобразование позволяет избавиться от избыточности непрерывного благодаря существованию функций, образующих полны® ортонормированные базисы пространства l2crd в случае, если ао= г и ьо= 1. Кроме того, ъ рамках дискретного преобразования при выполнении некоторых дополнительных требований удается построить быстрый алгоритм получения вейвлет-коэффициентов - так называемый анализ с переменным разрешением

С mul t i re*ol ut-i on anal ysi s5 Г 1 03 .

Анализ с переменным разрешением основывается на рассмотрении последовательности вложенных друг в друга аппроксимирующих функциональных пространств

V с V с V с V с V с . - - . С 4D

г 1 о - i - г

таких ч что

и V = Лю, V. - <0>. С5.1

jeZ J je2 J

Дополнительным требованием, непосредственно связанным с; понятием переменного разрешения, является то, что все пространства являются масштабированной версией центрального пространства т. е.

ГСхГ) е V «^»гсг'х3е V . со

} °

И наконец, требуется существование функции ф е vo> такой, что

множество <Фоп= фсх - пэ; п <=• г> образует ортонормированный базис в V . Тогда <Ф с>о = гГ*'гфсг~’х. - го, п е г> образует

О ^ ;,1Ч Г

ортонормированный базис для V , при любых .* = г Функцию ф называют- масштабной функцией.

Таким образом, подпространство V яаля. ?ся подпространством функций, аппроксимирующих ь2сю с разрешением г5 и переход от пространства V к пространству V соответствует переходу к большим масштабам, т. е. огрублению функщи с естественной потерей некоторой информации, содержащейся в исходной функции.

Основная идея анализа с переменным разрешением состоит в построении алгоритма, использующего разность информации, содержащейся в различных масштабах. Для любого j е г, определяется подпространство - ортогональное дополнение V в у ^.-

V = V # « , V х V и V» л. V , при \ * 1' ■ с 7 3

^ * J ) } Л Л 3 Г

Кроме того, из свойства сео следует, ЧТО 1/^1»= Ф^ , где все и -ортогональны друг другу. Очевидно, подпространствам также присуще масштабное свойство с Для любой последовательности замкнутых подпространств, удовлетворяющих условиям с 4-ео, существует вейвлет-функция г. такая, что множество <у = г" *у'Аусг~ -’х - к^, к е г> образует ортонормированный базис подпространства у. для любого _ъ Полное множество <> к; .).!< е г> есть ортонормированный базис для ь2сю.

Дополнительное упрощение, введенное Малла поз позволяет вычислять коэффициенты разложения исследуемой функции, используя дискретные фильтры и д . Фильтр можно подучить из

вложенности подпространств <- V и ТОГО, ЧТО ф

ортонорм1фованный базис подпространства V .•

п<=г;

отсюда 1-^= с ф. Ортонормированность масштабных функций

налагает на условие; Е ьГ(ьп+2к = 6ко с ^-символ Кронекера^.

Т>Е2

Другой фильтр зп связан с вейвлет-функцией ч> и вытекает из вложенности подпространств с V' = Е 9п^_1п- причем

существует простая СВЯЗЬ между Ь И 9Г1: дт =С-Оп или

иногда дг = С -15п Ь_ )+1+гм с подходящим выбором N е 2.

Анализ с переменным разрешением приводит к быстрой иерархической схеме для вычисления вейвлет-коэффициентов

исследуемой функіда г поз.

Начиная с коэффициентов разложения исследуемого сигналя для самого мелкого масштаба с самого высокого разрешениям .леї'ко ЬОДУЧЮЪ. используя СЙО , коэффициенты Сі * 1,2,3 .3, 110

формулам ООО И вейвлет....коэффициенты С ^ 0,1,2. АЛГОрИГМ

можно рассматривать как вычисление последовательно огрубляемых аппроксимаций г, определяемых ко&ффициентами с£, с одновременным учетом информации, теряемой при переходе к боле** котлов масштабу. Эта информация и отражается в вейвлет-копффщиентах а1 Ка практике наименьший возможный масштаб огцюделяе-'-оя числом дискретных значении і анализируемого» сигнала ; .удобно, чтобы х * г*', з е- При переходе от данного масштаба к следующему число вейвлет-коэффицдантов уменьшается в два рааа и процесс останавливается после конечного числа уровней. Так как масштабная функция хорошо локализована в прострзистьо. тс в качестве коэффициентов г" = <г,ф у можно взять сами доскре' ные значения исходной функции.

В терминах функциональных пространств этот алгоритм

соответствует разложению пространства V на сумму ортот оаальных подпространств вейвлет-функций последовательно огрубляемых масштабов; V ■- я ф н &.. . & ш фу .

О О і J-± .9 J

Хотя условия с 4-оо кажутся весьма сложными, г;р<е,лложе.ч' достаточно большое число функций. удовлелъеряюшкх требованиям аи&чиза с переменным разрешением, например, вейвлеты Мейера, имеющие оіраничен.чьт спектр в фурье-пространстве, зейвлегы

Баттл-Лемари. основанные на сплайн-функциях, и вейвлеты Цобеши

[61 - гладкие функции ОТЛИЧНЫ© ОТ нуля Лі'ВіЬ НЭ КОЖИНОЙ

интервале, поэтому соответствующе им фильтры Ь и -.} КШЮТ конечное число ненулевых членов. Приводам значения коэффициентов фильтров для функции Каара: ь = ь = а'1х2, и для Функции Лобегаи з-го порядка: ь = о.ээгеуовз, ь = о. всевеїш, ь = о. 45997750,

Ь = -О. 13Р01102, Ь => “0.08544127, Ь = С.С35г2&2б.

Я 4 5

БЕйВЛЕТ-ПАКЕТЫ. Описанный выше алгоритм оставляет открытым вопрос выбора самого базиса, ответ на который обычно не однозначен и зависит от характера исследуемого сигнала. При анализе аэрометрической информации весьма перспективным может оказаться использование предложенного недавно Викерхаузерсм из.1 метода войвлет- -пакетов. Пакеты расширяют класс базисных функций и предоставляют целые "библотеки" базисов, которые подстраивают-Г'-1 под конкретные масштабные и частотные особенности анализируе-„го сигнала. Грубо говоря, вейвлет-пакет есть множество функций из пространства ь2ср?:> с модулированными колебаниями, хорошо .локализованных как в пространстве, так и по частоте. С каждой функцией можно связать три параметра: масштаб, местоположение и частоту. Семейство вейвлет-пакетов получается посредством растяжения, переноса и мг-тул* ’материнского" вейвлета.

Построение г)учетов можно осуществить с использованием тех же дискретных фильтров ь и д, с помощью которых выполняется, анализ с переменным разрешением. Определим последовательность функций:

$ С

2г>

1

С гх - .р

си:

Ф СхМ

2Г1+ 1

С 2х

сз;

ФоС>:М МОЖНО отождествить С МЗСШТабНОЙ функцией ф. З Ї Сх) - с вейвлет-функцией у/. Определим функциональное пространство

Тогда множество функций §пСх-.р, j е z являем :,.я ортонормирован-ным базисом пространства с>п Доказано с i зз, что функции s cx-j>, n е z, n > о образуют ортонормированный базис l2cr_v Далее, введем пространство 6ko^, получаемое масштабным преобразованием о :

п

<5*0 = <f: f =■ Г:,,гк'2$ С£кХ ~ р, со е R> .

Г> 1 J n J

Множество функций сгкХ2$г1сг!х - p. j е z> образует вейвлет-пакет. Из этого пакета функций можно выделить большое число различных ортонормированных базисов пространства i Лк. Например, если взять в качестве' исходной функцию Хаара, то вейвлет-пакет будут образовывать функции Уолша различного порядка. Вейвлет-базис является частным случаем вейвлет-пакет базисов. Соответствующая ему последовательность подпространств

ИМееТ ВИД:

О Ф О Ф 60 Ф. . . Ф 6 О Ф. . . .

Oil 1

Пусть f = гсх:. - функция из пространства l2cr:> и . р «= z> коэффициенты проекции 1' на пространство самого мелкого

масштаба.- г = ci\2LX2$ са1'х - рох Тогда вводя обозначения

Р L, *-*

l'nk=C f , 2* Х2ф С2 х-рЮ , р е Z, О < )с :? L, О < n < S1' У , ПОЛУЧЭбМ

р п ^

реккурентные соотношения;

f 2„.k-l = у h f г.к J

р L j-zp j

С 14.'»

Используя эти соотношения легко получить все, возможные коэффициенты г^к, р е- г - соответствует положению базисной функции.

О 5- к •: ! , - масштабу с отметим, что здесь большим к соотвеї ствуг.т малые масштабы:-, о * п <- 2*''к - частоте. Коэффициеты разложения удобно представить в виде векторов <г'\ р е- ху, организованных в бинарное дерево. Корню дерева соответствует пространство а потомкам і.-к -го поколения - подшростра;* .-тва ••> - о. ... ,

Алгоритм вычисления векторов коэффициентов можно рассматривать как построение дерева начиная от корня: коэффициенты следующего уровня получаются из предыдущег; путем применения соотношений сізз и с 143. Ниже приведена схема этой процедуры для восьмиточечного сигнала: одинарная вертикальная

линия обозначает применение фильтра ь с формула і зз, а двойная вертикальная линия - фильтра д Ч4з. Для примера выделены двойной рамкой коэффициенты, соответствующие вейвлет-базису.

Из данного вейвлет-пакета можно выделить очень большое число оргснормированных вейвлет-базисов. Нужно из этого множества базисов выделить наиболее подходящий для разложения данного сигнала. Например, Викерхаузер г гзэ предложил быстрый алгоритм поиска базиса, минимизирующего некоторую ценовую функцию коэффициентов, например информационную энтропию. С практической точки зрения такой базис минимизирует число существенных коэффициентов ст.е. коэффициентов, которые превосходят некоторое пороговое значение:). Базисные функции, соответствующие этим коэффициентам, содержат почти всю информацию о сигнале и могут позволить выделить его характерные моды.

Первый этап этого алгоритма состоит в вычислении ценовой

функции хая всех узлов дерева. Затем, проходя дерево снизу вверх, отбрасываем нижележащие ветви, если значение ценовой функции отца меньше суммы ценовых функций сыновей Тогда для данного сигнала и ценовой функции получается вполне определенный ортонормировэнный базис.

Для сигнала, представленного в виде і - гк дискретных значений, алгоритм *икерхаузера ищет оптимальный базис из множества, содержащего г1 различных ортонормированных базисов. Однако этот алгоритм определяет далеко не все возможные базисы, но эффективный способ просматривания полного множества неизвестен. Кроме того, ‘оптимальный базис с точки зрения минимума энтропии не всегда является подходящим для выделения конкретных особенностей сигнала. Например, хорошо идентифицируя главные гармоники сигнала . такой базис может "игнорировать" мелкомасштабную структуру сигнала. Эту проблему можно решать введением некоторых эмпирических требований на вид базиса, например ограничивать максимальный используемый масштаб, или варьировать длину анализируемого сегмента сигнала.

НЕКОТОРЫЕ ПРИМЕРЫ. Объем статьи позволяет продемонстрировать возможности вейвлет анализа .лишь на нескольких характерных примерах. Первый пример не имеет отношения к аэрометрическим данным, но хорошо демонстрирует возможности подхода срис.г:.. Сигнал представляет собой сумму двух периодических составляющих, каждая с линейно меняющейся частотой. В нижней части рисунка приведен сам сигнал, верхняя часть - фазовая плоскость время-частота. Отметим, что наглядное представление результатов вейвлет преобразования является самостоятельной проблемой. В нашем случае при построении рисунков фазовая плоскость разбивается на прямоугольники <т = соответствующие

областям локализации базисных функций, В силу принципа неопределенности прямоугольники имеют равную площадь, ширина прямоугольника характеризует временной интервал, а высота -диапазон частот, интенсивность окраски пропорциональна квадрату коэффициента )^к, индексы р и п - определяют координаты центра прямоугольника, а к - его ширину. Ортогональность и полнота

базиса гарантируют, что фазовая плоскость будет покрыта прямоугольниками без перекрытий и пропусков. Базис для данного сигнала выбран та условия минимума энтропии с применением алгоритма быстрого поиска.

Второй пример ближе к проблемам экологии и соответствует сигналу, включающему две гармонические компоненты и три случайных выброса срис,зэ. Причем, если два из них хорошо видны на графике самой функции, то третий, имеющий такую же амплитуду, почти не заметен. В поле вейвлет-представления постоянной частоте соответствует горизонтальная периодическая структура, а локальному выбросу - вертикальная. К сожалению, качество печати не позволяет передать видимую на экране дисплея мелкомасштабную структуру рисунка.

Третий пример касается суточной изменчивости концентраций примеси в атмосфере. Известно п. аз, что для окиси углерода почти во все сезоны наблюдаются два вида суточной изменчивости» либо с одним, либо с двумя максимумами. Сигнал на рис.4 сформирован из фрагментов» точно соответсующих этим двум типам изменчивости, взятым из си, плавно сменяющих друг друга. Задача состоит в разделении интервал-"' времени, в которые доминирует тот или иной тип зависимостей.

Приведенные примеры лишь иллюстрируют возможности подхода. Остается много проблем как общего плана, связанных с выбором базиса, оптимизацией алгоритмов, представлением результатов, так и специфических, возникающих при работе с аэрометрическими данными. К последним относится в первую очередь проблема пробелов в выборках - замеры на постах наблюдения проводятся, как правило, три раза в сутки, в 7, 13 и 19 часов и не проводится в ночное время, а также в выходные и праздничные дни. В то же время именно этот аппарат дает надежду на выделение из достаточно сложных временных последовательностей данных Это могут быть как отдельные метеопараметры или концентрации примесей, так и интегральные показатели состояния атмосферы характерных временных структур с целью их использования для дальнейшего прогноза состояния системы.

Рис. 2

^/VVVYVVVVVVЛ/^\/VЛ/^;

Рис, 3

J]xЛAйf^fЩAЛ^A^AAAAAA^J\AA^f^MA/Jl^AA

Рис. 4

Литература

і . Безуглая Э.Ю. Мониторинг состояния загрязнения атмосферы. Л.: Гидрометеоиздат, 1986.

Берлянд М.Е. Прогноз и регулирование загрязнений атмосферы. Л.: Гидрометеоиздат, Г985.

3. Зимин В.Д. Иерархическая модель турбулентности // Изв. АН ССОР. Физика атмосферы и океана. 1981. N-12. С.1265-1273.

4. Зимин В.Д., Фрик П.Г. Турбулентная конвекция. М.:Наука, 1988.

Применение В е й е? л е т - а н а л и з а

з. Фрик П. Г. Вей&дэт-анализ и иерархически© модели

турбулентности, - ИМСС УрО РАН. Пермь» 1392.

6. Daubechies I. Tan lectures on wavelets. .// CBHS Lecture Nr у tes IT г i es , SI AM. 19Q1 .

7. Farge M. Wavelet transform and their applications to

turbulence. //Annu. Re-» Fluid Mech. X©92. N.£4 P.3Q5-4S7.

9. Frick P. G. , Zimin V. D. Trans. Int.Conf. "Wavelets,

fractals and Fourier transforms: new developments and new

applications", Cambridge, Dec. 1S90. - Oxford Press, 1991.

"5. Gr ossmann A. , Morlet J. Decomposition of Hardy functions into square integrable wavelets of constant shape. // SIAM J. Math. Anal. 1984 V. IS. N. 4 P. 723-736.

10. Mall at S. A theory for multiresolution signal

decomposition: the wavelet representation. // IEEE, Trans. on

Pattern Anal. Machine Intell. 1989. N 2. P. 7.

11. Meneveau C. Analysis of turbulence in th<? orthor.oriral wavelet representation. // J. Fluid Mech. 1991.

12. Meyer У. Ondeiettes et operateurs. Paris: Hermann. 1990.

13. Wi ckerhauser М. V. INRIA 1 ectures on wavelet packet

algorithms. У У Roquencourt, Jun. 17-21. 1991. P. 31-99.

Институт механики СПЛОШНЫХ сред