УДК 519.6
МНОГОКОМПОНЕНТНАЯ МОДЕЛЬ ВРЕМЕННОГО РЯДА
О.В. Мандрикова (КамчатГТУ)
Предложен новый подход к построению моделей временных рядов со сложной структурой, в основе которого лежит процедура разложения сигнала по базисным функциям. В качестве базисных функций рассматриваются ортогональные вейвлеты. Рассмотренная в работе конструкция позволяет выявить и исследовать те особенности структуры данных, которые не попадают в область традиционных методов анализа временных рядов.
The author offers new approach to composing temporal series models with complex composition. The basis of this approach is the procedure of splitting the signal into basic functions -orthogonal wavelets. The construction offered by the author allows to trace those features of data structure that are not covered by traditional methods of temporal series analysis.
Введение
В большинстве случаев, когда требуется выделить регулярные компоненты сигнала и построить прогноз, удобным является метод идентификации модели временного ряда на основе класса моделей авторегрессии проинтегрированного скользящего среднего (АРПСС) [1]. Данный метод чрезвычайно популярен во многих приложениях. Практика подтвердила его мощность и гибкость, но он имеет ограничения как на возможность его использования для отдельных временных рядов, так и на выявляемые при этом закономерности. По этой причине данные методы не могут быть использованы для сигналов со сложной структурой (например, сигналов речи и изображений, сложных природных сигналов). Получившие бурное развитие в настоящее время детерминированные методы аппроксимации, базирующиеся на разложении функции по базисам, объединили теорию аппроксимации и обработку сигналов [2]. Разработанные на их основе адаптивные алгоритмы дают возможность на основе небольшого числа параметров получить аппроксимации требуемой точности и, что самое главное, могут быть использованы для широкого спектра сигналов с различной структурой. Во многих прикладных задачах требуется идентифицировать модель сигнала с целью предсказания его будущих значений. Идентификации, как правило, подлежат сигналы, включающие в себя различные по длительности и форме локальные особенности, которые содержат полезную информацию и не могут быть отфильтрованы как шум. Примерами могут служит сигналы регистрации геофизических и геохимических параметров, содержащие особенности в виде всплесков, серии пиков различной амплитуды и частоты, ступенеобразные особенности. Последние делают невозможным непосредственное применение к ним традиционных методов анализа временных рядов.
Многокомпонентная модель временного ряда
С целью расширения возможностей применения математического аппарата к обработке статистических данных введем в рассмотрение следующую математическую конструкцию: измеряемую величину f определим как
f О) = Xlfl (t) + Х2 f2(t) + ... + Xnfn (t) , (1)
где x1, x2,..., xn - параметры; t - время.
I. На функции f наложим выполнение условия ортогональности относительно величин t. с весовыми коэффициентами gt = 1/о2 (условие 1):
Z gf (f.) f (f.)=5*.
i
Введем предположение, что ошибки измерений являются белым шумом (условие 2). При выполнении условий 1, 2 ковариационная матрица для коэффициентов x1, x2,..., xn является единичной матрицей, т. е. все эти коэффициенты не коррелируют между собой.
II. С целью применения данной конструкции для сигналов с различной структурой будем предполагать, что функции / могут быть различными, в частности иметь разную форму (условие 3). Компоненты / (I) идентифицируем как / (I) = I - кфЛ (I), / е Ь2(К, где функции фЛ -
к
г 2 /
базисные функции пространства Ь (К) (условие 4).
Математическую конструкцию (1) с учетом указанных условий 1-4 назовем многокомпонентной моделью временного ряда (ММВР).
Введенная математическая конструкция раскрывает широкие возможности при решении различных задач аппроксимации данных. Используя соотношение (1) с учетом условий 1 и 2, мы рассматриваем сигнал как реализацию случайного вектора. В свою очередь, условие 4 включает в себя детерминированный подход, который базируется на разложении функции по базисным функциям. Поскольку современные детерминированные методы аппроксимации функции предоставляют широкий спектр базисных функций [2], это даст возможность идентифицировать те особенности структуры данных, которые не могут быть описаны на основе традиционных стохастических методов. Компоненты модели имеют более простую структуру, чем исходный сигнал. Оценка параметров с,к его компонент в зависимости от их структуры и имеющейся априорной информации может быть выполнена на основе традиционных стохастических методов.
Способ идентификации компонент ММВР
В качестве метода идентификации компонент ММВР будем использовать вейвлет-преобразование. Любой сигнал конечной энергии / может быть разложен по ортогональному
базису пространства Ь (К) . Вейвлет ¥, п является единственной базисной функцией этого пространства, позволяющей восстановить дискретный сигнал на основе неполной выборки. Вейвлет-конструкция обладает важным преимуществом - она включает в себя широкий спектр базисных функций различной формы. Это свойство обеспечивает выполнение условия 3 ММВР и дает возможность использовать эту конструкцию для большого числа сигналов с различной структурой.
На основе вейвлет-конструкции сигнал представляется единственным образом в виде суммы разномасштабных составляющих:
V/ е Ь2(Я)3!/(I) = ... + уч(|) + у0(I) + ) + ..., у е * е % , (2)
где у - составляющая / . Каждая такая составляющая дает локализованную частотновременную информацию об / в определенном частотном диапазоне.
С учетом этого свойства определим функции / в конструкции (1) следующим образом:
/ (') = 1 с, ,, к ¥, к (I), / е Ь2 (К) . Каждая / е Ь (К) может быть записана как
, ,к
/«) = I с,к * (I) = 11М1 , * „ (I), (3)
где
Эти бесконечные ряды, называемые вейвлет-рядами, сходятся в Ь2(К) [3, 4]. Конструкция (3)-(4) определяет дискретное вейвлет-преобразование (ДВП).
Утверждение (о допустимости использования конструкции ДВП в качестве метода идентификации компонент ММВР). Конструкция ДВП осуществляет разложение функции / еЬ2(К) на ортогональные компоненты /(I) = х /1(1) + х2/2(1) +... + хп/п(I) +... и обеспечивает
1 I 1 I
выполнение условия ММВР (ґі)/к (ґі) = 5,к при условии, что X. = ^ с2к (і,) , где
'■ С< I , к , ,
с. к - коэффициенты вейвлет-разложения.
Доказательство. Структура разложения 1}(К), порожденная вейвлетом ¥е£2(^), имеет вид
И(К) = Хж : = ... +Ж_1 + Ж° + Ж + ...,
(5)
.ех
где Ж : = clos 2 (т ; к е X).
^ . Ь2(К)\ . ,к’ /
Функция / при этом представляется в виде суммы компонент (см. (2)):
V/ е Ь2(Я)3!/(I) =... + у_,(0 + ^(0 + vl(l) +...,V,. е Ж,г е X. (6)
Каждая компонента V. имеет единственное представление в виде вейвлет-ряда (см. (5)):
V. = Х.т и (I). (7)
кеХ
Если Т - ортогональный вейвлет, то подпространства Ж. в соотношении (7) - взаимно ортогональны (см. (6)), т. е.
(у., V.) = ^ (8)
где V. е Ж , V. е Ж .
^ ’ г г
V.
Определим в ММВР функции / как / = —. Тогда получим:
' ' х.
- если г * I, то
- если г = I, то
Хз- / (I.) / (I.) = °;
/ °2
Хтт/ (I.) / (I.) = Етг /ЧО) = Х^Т X с2к (I. ) =
(9)
(1°)
где коэффициенты с. к определяются из соотношения (7) V . = Х с. ,к¥. ,к (I).
кеХ
г , У/2
Положив х =
Х4- X с2к (I.)
из соотношений (8)-(1°) следует выполнение условия
=«.. ммвр.
Таким образом, мы получили конструкцию (5)-(7) разложения функции / на компоненты, которые в силу соотношения (1°) ортогональные, и обеспечили выполнение условия ММВР
X 8/ (I.) / (I.)=5 *; ё = 1/°2.
г
Для обработки дискретных данных в вейвлет-теории используется конструкция кратномасштабного анализа (КМА). Данная конструкция имеет два важных преимущества:
1. Она позволяет выполнять обработку дискретных данных. В этом случае, имея дискретные значения функции / (т. е. значения функции на сетке с разрешением 2-1), в качестве пространства выборки рассматривают подпространство V. = closL2^R (ф(211 - к)) , где функция ф ,
порождающая подпространства V. своими сдвигами и растяжениями, называется масштабирующей функцией. Она играет особую роль в конструкции КМА. От свойств этой функции зависят такие важные свойства вейвлета Т , как величина носителя и гладкость [2-4].
2. На конструкции КМА базируется алгоритм быстрого вейвлет-преобразования сигналов конечной длины, который требует только 0(N операций для сигналов длиной N .
Идентификация компонентов ММВР на основе конструкции КМА
Предположим, что сигнал / известен с разрешением 1 (т. е. у = 0). Тогда базисным пространством для него будет У0 = є/о,? 2 (ф(20 ґ _ к)), к є 2 . Выполним разложение сигнала /0 на
Ь (Л )
две компоненты: /0(ґ) = g_1(t) + /_1(ґ) . Ортогональная проекция исходного сигнала /0(ґ) на подпространство У_1 даст нам компоненту /_1(ґ). Она задается набором скалярных произведений /0(ґ) с функциями из ортобазиса У_1, т. е. величинами
= (/0(0,2_1/2ф(-2_к)^ = ^/пф(ґ_и),£Ьф(ґ_= £Ь/и+,. (11)
В качестве компоненты g_1(t) рассматривают компоненту сигнала, ортогональную к пространству У1. Получают схему разложения У0 = У1 ® Ж1. Ортобазисом пространства Ж1 явля-
=л £ цкф(ґ _ к), коэффициенты цк
к
имеют вид дк = (_Х)кН1_к [2]. Проекция сигнала на подпространство Ж_х задается набором скалярных произведений /0 (ґ) с функциями ^ :
что равносильно свертке с фильтром q* и прореживанием вдвое.
Описанная вычислительная процедура называется алгоритмом разложения. Она работает на любом масштабе и позволяет выполнить разложение сигнала на компоненты по схеме V. = V.-1 © Ж.-1. Если эту процедуру разложения выполнить для сигнала /й(!) т раз, то получим следующую схему разложения:
V = Ж. © Ж, ©... © Ж © V , (13)
° -1 -2 -т -т ’ V /
где Ж., V. - пространства с разрешением 2-].
В результате исходный сигнал /° будет иметь единственное представление в виде конечной суммы компонент:
Л 0) = 8-1 0) + ё-2</) + ... + ё-т 0) + /-т (I). (14)
Каждая компонента единственным образом определяется последовательностями коэффициентов С и й . :
' /.(I) = Х с.к ф( 2.I - к )еУ,,
к
С := {.} к е Х,
' 8. (I) = Х . Т( VI - к )еЖ,
<_ к (15)
:={.}, к е X.
Коэффициенты й. являются детализирующими компонентами сигнала. Коэффициенты С соответствуют сглаженной компоненте сигнала.
Коэффициенты С позволяют выделить полезную информацию о сигнале, находящуюся в низкочастотной его составляющей. Постепенное уменьшение масштабного параметра . позволяет фокусироваться на локальных структурах сложного нестационарного сигнала и исследовать его структуру. Поэтому коэффициенты й. обеспечивают идентификацию тех особенно-
ется набор вейвлет-функций 2 1/2 -2- _ т|1 , где 2
стей структуры данных, которые не могут быть описаны на основе традиционных параметрических методов моделирования, но содержат полезную информацию о процессе в широком спектре приложений.
Аппроксимация компонент ММВР на основе вейвлетов
Приближение / М векторами, индексы которых принадлежат множеству ІМ , имеет вид
/М = £( Л 'т)'т .
тєІт
Погрешность аппроксимации - это сумма квадратов модулей скалярных произведений с векторами, индексы которых не принадлежат множеству ІМ :
АМ ]=|| / _ /м\ I' = Ц< /. О I’.
||2
тИт
Чтобы минимизировать погрешность, следует выбирать 1М таким, что М векторов с индексами из 1М имеют наибольшие модули скалярного произведения.
Очевидно, амплитуды скалярного произведения |(/, ет^| связаны со структурой сигнала. Используя в качестве базисных функций вейвлеты ТЛк и выбирая апостериори т векторов Т3,к, зависящих от свойств функции / , мы получаем приближение / т векторами, индексы которых принадлежат I (см. (13), (14)):
/(I) =Е ^ ,кТ.,к + /-т (I). (16)
. ,ке1
В вейвлет-теории доказано, что амплитуды скалярного произведения |</■Т ">| локально
возрастают там, где сигнал негладкий [2-4]. Поэтому нелинейная аппроксимация, которая содержит наибольшие скалярные произведения, эквивалентна построению адаптивной приближающей решетки.
Аппроксимация сглаженной компоненты сигнала на основе модели АРПСС
Решая задачу построения модели сигнала на основе конструкции ММВР и используя в качестве метода идентификации ее компонент КМА, мы получаем помимо мелкомасштабных компонент 8-1 (I) сглаженную (низкочастотную) составляющую сигнала /-т (I) (см. (14)). В идеале хочется, чтобы низкочастотная компонента /-т (I) позволяла получить наилучшее приближение функции / , а с помощью высокочастотных составляющих ё-1 (I) идентифицировались лишь отдельные локальные высокочастотные ее характеристики. Поэтому наилучшую аппроксимацию высокочастотных составляющих для большинства сигналов может обеспечить нелинейная вейвлет-аппроксимация, предложенная в соотношении (16). Поскольку сглаженная компонента разложения /т (I) имеет более простую структуру, чем сам сигнал, то это в большинстве случаев позволяет использовать для оценки ее параметров традиционные параметрические методы моделирования временных рядов. Например, в случае использования модели авторегрессии модель сигнала будет иметь вид (см. (16)):
/(I) = £ й.,кТ.,к +£ у,, (17)
. ,ке1 /=1
где = Vлс,к; у1,..., ур - коэффициенты авторегрессии; V - оператор взятия разности; й - порядок разности; Тк (I) = 21/2Т(2Ч - к) - базисный вейвлет; Т3,к - соответствующий ему двойственный вейвлет.
Сглаженная компонента сигнала /-т (I) аппроксимируется АР-моделью. Для идентификации АР-модели используют итеративный подход [1, 5], а именно:
- на основе грубых методов идентификации получают грубые предварительные оценки параметров модели;
- далее используют грубые оценки как начальные значения в более точных итеративных методах оценивания параметров;
- диагностические проверки позволяют выявить возможные дефекты подгонки и диагностировать их причины. Если такие дефекты не выявлены, модель готова к использованию. В противном случае итеративные циклы идентификации, диагностические проверки повторяются до тех пор, пока не найдено подходящее представление модели.
В нашем случае перечисленные этапы идентификации АР-модели образуют промежуточный этап общей последовательности этапов идентификации модели (17). С учетом сложной структуры аппроксимируемых данных при подгонке модели (17) этапы идентификации АР-модели должны быть выполнены для нескольких сглаженных компонент сигнала, соответствующих разным уровням вейвлет-разложения. На основе результатов диагностических проверок полученных компонент модели и с учетом целей исследования должна быть идентифицирована окончательная модель. Диагностическая проверка пробной модели включает в себя:
- оценку погрешности компоненты f(t) модели;
- оценку погрешности компоненты f2(t) модели.
Литература
1. Никифоров И.В. Последовательное обнаружение изменения свойств временных рядов. -М.: Наука , 1983.
2. StephaneMallat. A Wavelet tour of signal processing / Пер. с англ. - М.: Мир, 2005.
3. Ingrid Daubechies. Ten Lectures on Wavelets / Пер. с англ. - Ижевск: Регулярная и хаотическая динамика, 2001.
4. Charles KChui. An Introduction to Wavelets / Пер. с англ. - М.: Мир, 2001.
5. Бокс Дж., Дженкинс Г. Анализ временных рядов: прогноз и управление / Пер. с англ. -М.: Мир, 1974.
6. Гидрогеохимические предвестники землетрясений / Под ред. Г.М. Варшала. - М.: Наука, 1985.
7. Беляев А.А. Особенности радоновых прогнозных признаков землетрясений // Геохимия. -2001. - № 12. - С. 1355-1360.
8. Мандрикова О.В. Непараметрическое моделирование сейсмических данных с использованием сплайн-вейвлетов // Изв. СПбГЭТУ. Сер. «Информатика, управление и компьютерные технологии». - 2002. - № 4. - С. 77-81.
9. Фирстов П.П., Филиппов Ю.А., Мандрикова О.В. Предвестниковые аномалии сильных землетрясений в динамике подпочвенного радона на Петропавловск-Камчатском геодинамиче-ском полигоне в 1997-2001 гг. // ДАН. - 2003. - Т. 389. - № 6. - С. 810-813.