ISSN 0236-235X (P)
ISSN 2311-2735 (E)

Публикационная активность

(сведения по итогам 2017 г.)
2-летний импакт-фактор РИНЦ: 0,500
2-летний импакт-фактор РИНЦ без самоцитирования: 0,405
Двухлетний импакт-фактор РИНЦ с учетом цитирования из всех
источников: 0,817
5-летний импакт-фактор РИНЦ: 0,319
5-летний импакт-фактор РИНЦ без самоцитирования: 0,264
Суммарное число цитирований журнала в РИНЦ: 6012
Пятилетний индекс Херфиндаля по цитирующим журналам: 404
Индекс Херфиндаля по организациям авторов: 338
Десятилетний индекс Хирша: 17
Место в общем рейтинге SCIENCE INDEX за 2017 год: 527
Место в рейтинге SCIENCE INDEX за 2017 год по тематике "Автоматика. Вычислительная техника": 16

Больше данных по публикационной активности нашего журнале за 2008-2017 гг. на сайте РИНЦ

Вход


Забыли пароль? / Регистрация

Добавить в закладки

Следующий номер на сайте

4
Ожидается:
16 Декабря 2018

Алгоритм выделения тренда зашумленных больших временных рядов

Algorithm for identify trends for large noisy time series
Статья опубликована в выпуске журнала № 4 за 2011 год. [ на стр. 36 – 40 ][ 13.12.2011 ]
Аннотация:В работе предложен эффективный и точный метод нахождения тренда, фильтрации шумов и сглаживания боль-ших нестационарных временных рядов. Показано, что тренды таких рядов могут быть определены при помощи вейвлет-пакетного преобразования совместно с сингулярно-спектральным анализом.
Abstract:In this work efficient and accurate method for identify trends, denoising and smoothing for large nonstationary time series are provided. It is shown that the trends such large series can be determined the wavelet – packet transformation in combination with singular-spectrum analysis.
Авторы: Бахмутский М.Л. (mbakhmut@mail.ru; ) - НИИСИ РАН, г. Москва, , , кандидат физико-математических наук
Ключевые слова: сингулярно-спектральный анализ, вейвлет, временные ряды, выделение тренда, алгоритм
Keywords: singular-spectrum analysis, wavelet, time series, identify trend, algorithm
Количество просмотров: 13104
Версия для печати
Выпуск в формате PDF (5.83Мб)
Скачать обложку в формате PDF (1.28Мб)

Размер шрифта:       Шрифт:

Выделение тренда и сглаживание больших, порядка десятков тысяч значений, временных рядов является актуальной задачей. Для ее решения , предлагается применять сингулярный спектральный анализ в сочетании с предварительным применением к исходному временному ряду вейвлет-пакетного разложения.

Сингулярный спектральный анализ

При изучении временных рядов, помимо анализа Фурье и вейвлет-анализа, широко применяется метод сингулярно-спектрального анализа (ССА) [1]. Данный метод основан на построении по исходному временному ряду, по так называемой траекторной матрице, ее аппроксимации при помощи сингулярного разложения матриц. Поэтому его можно рассматривать как исследование временного ряда хорошо известным в статистике методом главных компонент (другое название – метод естественных ортогональных функций). Для построения траекторной матрицы временной ряд  последовательно, со сдвигом на один отсчет, проходится окном размером m отсчетов и строится множество векторов задержки размерности m  [2], которые составляют траекторную матрицу размерности m*n

A=,              (1)

где  – значения ряда в моменты времени t=1, 2 ,…, N, n=N-m+1, m

Для траекторной матрицы находится ее сингулярное разложение, то есть матрица А представляется в виде A=U, где U – матрица размерности m*n с ортонормальными столбцами, U=; V= – матрица размерности n´n с ортонормальными столбцами; =diag () – матрица размерности m*n с неотрицательными диагональными элементами (сингулярными числами). Матрица аппроксимируется при помощи r троек , отвечающих r наибольшим сингулярным числам. Иными словами, матрица A заменяется матрицей :

=,                 (2)

где = – внешнее произведение i-го столбца матрицы U и соответствующего столбца матрицы V;  – простая матрица ранга 1 с единичной нормой. Временной ряд восстанавливается путем диагонального усреднения матрицы , то есть матрицы, в которой оставлены только r более гладких компонент.

=                  (3)

Если параметр m выбран таким, что m, можно решить спектральную задачу для матрицы C=A*=U*(*)*.

Матрица U – это матрица собственных векторов матрицы C, а сингулярные числа матрицы A являются положительными корнями собственных чисел матрицы C. Умножив скалярно матрицу  на A, получим матрицу главных компонент: Y=*A=, где , i=1, 2 ,…, m – строки длины n.

Подпись:  Рис. 1. Исходный временной ряд 130 816 значений давления, измеренных с интервалом 1 минута Рис. 2. Исходный ряд (черная линия) и его сглаживание (светлая линия) методом ССА с движущимся окном в 100 отсчетов Рис. 3. Исходный ряд и его сглаживание методом ССА с движущимся окном в 1 000 отсчетов Рис. 4. Вторичное сглаживание исходного ряда методом ССА (сплошная черная линия – окно 100 отсчетов, сплошная светлая – окно 500 отсчетов, светлые точки – 5 000 отсчетов)Тогда = и временной ряд восстанавливается по формуле (3). Этот подход довольно широко используется в задачах нахождения периодических зависимостей во временных рядах наблюдений, фильтрации шумов и сглаживания временных рядов, так как использование сингулярного (спектрального) разложения матрицы позволяет выделить наиболее значимые составляющие ряда и отсеять случайные возмущения. Однако при непосредственном применении этого подхода к обработке очень больших рядов наблюдений возникают трудности. В изложенном методе ССА есть два существенных параметра – размерность окна, то есть размерность траекторной матрицы, и количество сингулярных троек (главных компонент), используемых для аппроксимации траекторной матрицы. Количество используемых главных компонент в данной задаче можно определять интерактивно.

Выбор длины скользящего окна определяется основной идеей метода ССА, то есть соображением о том, что, если исходный временной ряд имел какую-то структуру, то и его отрезки (столбцы траекторной матрицы) имеют такую же структуру. Если стоит задача выделения составляющих с известным периодом, размер окна равен или кратен этому периоду. Поэтому в задаче о выделении глобального тренда и сглаживании ряда рекомендуемый размер окна примерно равен половине длины ряда. При этом выбор размера движущегося окна оказывается очень существенным. На рисунке 1 изображен временной ряд, полученный в ходе реального эксперимента по гидродинамическому прослушиванию скважин. На рисунках 2 и 3, иллюстрирующих влияние выбора размера окна, приведены результаты сглаживания (белая линия) исходного временного ряда (черная линия) с использованием окон в 100 и 1 000 отсчетов соответственно. Были проведены численные эксперименты по последовательному применению метода ССА. Исходный ряд сглаживался методом ССА с использованием окна размером 100, 500, 5 000 отсчетов. На рисунке 4 приведено сопоставление результатов последовательного сглаживания для различных окон (100, 500 и 5 000 отсчетов). Из этих расчетов следует, что для надежного и корректного нахождения тренда временного ряда необходимо использовать траекторную матрицу с длиной столбца, равной примерно половине длины ряда. Заметим, что для вычисления сингулярного разложения матрицы размерности порядка 5 0005 000 на компьютере c CPU Intel Core i5-750 (2,67 GHz) потребовалось несколько суток непрерывной работы. Это вполне понятно, так как вычисление SVD-разложения матрицы размерности m*n требует порядка 4*  операций [3]. Отсюда следует, что непосредственно применить метод ССА к выделению тренда и сглаживанию больших временных рядов не удается, так как невозможно вычислить сингулярное разложение матриц размерности порядка 65 000*65 000 элементов, во всяком случае, для компьютеров традиционной последовательной архитектуры. Поэтому предлагается сочетать метод ССА и применение вейвлет-пакетов.

 

Обобщение метода ССА на случай больших временных рядов

Допустим имеется большой, порядка десятков тысяч отсчетов, временной ряд длины N. Будем полагать без потери общности, что N=. Пусть N=2*M. Вместо временного ряда X= = рассмотрим ряд Z, состоящий из двух временных рядов S и D. Z=, где S= и D=.

Переход от исходного ряда X к ряду Z происходит следующим образом:

, i=1, 2, …, M.      (4)

Обратное преобразование от S и D к X имеет следующий вид:

 , i=1, 2, . .., M; (5)

компоненты ряда S представляют собой локально сглаженные, усредненные значения, а компоненты D – уточняющие коэффициенты для элементов ряда X.

Фактически на основе (4) раскладывается ряд X по вейвлетам Хаара [4].

Преобразования (4), (5) являются взаимно-од­нозначными, как легко убедиться непосредственно, и сохраняющими нормы. То есть, если , , , то .

Таким образом, пренебрежение какой-либо величиной в преобразованном временном ряде означает (после обратного преобразования) пренебрежение такой же малой величиной в исходном. На следующем этапе применяем введенное преобразование к рядам S и D. То есть от ряда Z перейдем к ряду  по тем же правилам, что и при переходе от ряда X к ряду Z

 

.

Обратный переход:

Подпись:  Рис. 5. Самая гладкая часть 8-го уровня вейвлет-преобразования Рис. 6. Следующий после самой гладкой части массив поправочных коэффициентов 8-го уровня вейвлет-преобразования Рис. 7. Следующий (2-й после самой гладкой части) массив поправочных коэффициентов 8-го уровня вейвлет-преобразования

 , здесь i=1, 2, …,, .

Все то, что справедливо при переходе от ряда X к ряду Z, справедливо и для перехода от ряда Z к ряду . Ряд SS является временным рядом локально сглаженных, усредненных значений элементов ряда S, а ряд DS состоит из уточняющих коэффициентов. Аналогичное справедливо и для рядов SD и DD. При этом ряд SS описывает самую гладкую часть исходного временного ряда. Сохраняя взаимно-однозначное соответствие, удалось перейти от исходного временного ряда к ряду, состоящему из двух временных рядов (гладкой части и поправок), а затем к ряду, состоящему из 4 временных рядов (гладкой части и поправок различной гладкости). Последовательно применяя этот прием (называемый вейвлет-пакетным разложением [4]) k раз получаем временной ряд, состоящий из временных рядов, имеющих длину . Для каждого из этих  временных рядов построим свою траекторную матрицу (1) размерностью m*n, где . Вычислим сингулярные разложения этих матриц. Заметим, что эта часть алгоритма может быть эффективно выполнена как на традиционных последовательных, так и на параллельных системах. Анализируя сингулярное разложение траекторной матрицы наиболее гладкой части, определяем количество сингулярных троек (или главных компонент), описывающих тренд, и тем самым определяем тот порог, ниже которого пренебрегаем сингулярными числами всех траекторных матриц. После этого формируем аппроксимации (2) траекторных матриц, по формулам (3) восстанавливаем временные ряды и, последовательно применяя (5), получаем сглаженную версию исходного временного ряда.

Обратимся к примеру. Исходный ряд состоит из 130 816 отсчетов, то есть N=130 816, k=8, =511. На 8-м уровне вейвлет-пакетного разложения имеется 256 временных рядов длиной 511 значений, каждый из которых соответствует (как гладкая часть и поправочные коэффициенты) всему исходному временному ряду. На рисунке 5 приведен график гладкой части вейвлет-разложе­ния, на рисунках 6 и 7 следующие 2 временных ряда поправочных коэффициентов. Заметим, что график на рисунке 5 внешне похож на график рисунка 4, соответствующий последовательному двойному сглаживанию методом ССА с окном в 100 отсчетов, поскольку последовательное применение первой формулы преобразования (4) можно рассматривать как своеобразное усреднение в окне 256 отсчетов.

Подпись:  Рис. 8. Первые 20 сингулярных чисел из графика рисунка 10 Рис. 9. Исходный ряд (черная линия) и тренд, найденный с использованием 1-й сингулярной тройки Рис. 10. Тренд первого ряда вейвлет-коэффициентов. Построен при использовании 10 членов разложения (2)Подпись:  Рис. 11. Тренд первого ряда вейвлет-коэффициентов. Построен при использовании 11 членов разложения (2). Видно появление шумовой компоненты Рис. 12. Исходный ряд (черная линия) и тренд, найденный с использованием 10 сингулярных троек (10 членов аппроксимации (2)) Рис. 13. Исходный ряд (черная линия) и тренд, найденный с использованием 50 сингулярных троек. Шумовая компонента усиливаетсяНа рисунках 6 и 7 приведены следующие после самой гладкой части (рис. 5) массивы вейвлет-ко­эффициентов. Легко заметить, что их нерегулярность увеличивается, то есть растет влияние шумовой компоненты. Однако просматривается также закономерность. Влияние шума на графике на рисунке 7 больше, чем на рисунке 6, и тем более, чем на рисунке 5. На рисунке 8 приведен график 20 наибольших значений сингулярных чисел для 1-й и 2-й траекторных матриц. Предположим, что в разложении (2) решено ограничиться только первым членом. Тогда, отбрасывая все члены в разложении (2) для всех рядов коэффициентов, сингулярные числа которых меньше использованного, и k раз проделав обратное вейвлет-преобраз­ование, получим тренд (светлая линия) на рисунке 9. Последовательно добавляя в разложение матрицы (2) новые члены, повышаем детальность воспроизведения временного ряда. На рисунке 10 приведен график первого ряда вейвлет-коэффици­ентов, восстановленного с использованием 10 первых членов в (2). А на рисунке 11, где приведен тот же график, но восстановленный с помощью 11 первых членов (2), ясно видна появляющаяся шумовая компонента. Таким образом, для выделения тренда исходного временного ряда восстанавливаем первый ряд вейвлет-коэффици­ентов при помощи 10 членов разложения (2), во всех других рядах вейвлет-коэффициентов в разложениях типа (2) оставляем только те члены, сингулярные числа которых не меньше 10-го сингулярного числа из SVD-разложения первого ряда вейвлет-коэффициентов. По формуле (3) восстанавливаем ряды и, проделав обратое вейвлет-преобразование, получаем сглаженный тренд (рис. 12). Рисунки 11 и 13 иллюстрируют дальнейшее повышение числа членов в (2), ведущее к повышению точности аппроксимации исходного временного ряда и ко все большему появлению шумовых компонент как следствию излишней точности аппроксимации.

Сравним стоимость SVD-разложения траекторной матрицы исходного ряда в предположении равенства длины движущегося окна половине длины ряда и  траекторных матриц предлагаемого алгоритма. Стоимость SVD-разложения матрицы размерности m*n примерно равна  операций [3]. Если длина движущегося окна примерно равна половине ряда, то  и стоимость SVD-разложения  операций. Для данного алгоритма стоимость SVD-разложения  траекторных матриц с примерно равна  операций. Таким образом, стоимость вычисления SVD-разложения при ис-

пользовании предлагаемого алгоритма уменьшается в  раза.

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

Автор благодарен С.Г. Вольпину, предоставившему экспериментальный материал.

Литература

1. Vautard R., Yiou P., Ghil M. / Physica D. 1992. Vol. 58.

2. Broomhead D.S., King G.P. / Physica D. 1986. Vol. 20.

3. Деммель Дж. Вычислительная линейная алгебра. М.: Мир, 2001.

4. Столниц Э., ДеРоуз Т., Салезин Д. Вейвлеты в компьютерной графике. Ижевск: НИЦ «Регулярная и хаотическая динамика», 2002.


Постоянный адрес статьи:
http://swsys.ru/index.php?page=article&id=2908
Версия для печати
Выпуск в формате PDF (5.83Мб)
Скачать обложку в формате PDF (1.28Мб)
Статья опубликована в выпуске журнала № 4 за 2011 год. [ на стр. 36 – 40 ]

Возможно, Вас заинтересуют следующие статьи схожих тематик: