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

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

При изучении временных рядов, помимо анализа Фурье и вейвлет-анализа, широко применяется метод сингулярно-спектрального анализа (ССА) [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.

Тогда = и временной ряд восстанавливается по формуле (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

.

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

 , здесь 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 отсчетов.

На рисунках 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.