Analysis and elimination of noise components in time series with variable pitch

Sklyar Alexander PhD in Technical Science Associate Professor, Department of Applied Mathematics, Russian University of Technology 119602, Russia, g. Moscow, pr-t Vernadskogo, 78 askliar@mail.ru Other publications by this author

1. Введение

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

При анализе данных временного ряда и прогнозировании на его основе возникает множество задач, связанных с выделением трендовой, колебательной и составляющих ^[1,2.3,4]. Для выявления таких рядов используются различные методы в зависимости от характера данных ^[5,6,7,8]. В то же время анализ и обработка зашумленных данных вызывает значительные трудности. Возникает задача устранения, по возможности, такого шума. Для устранения такого шума используются различные методы сглаживания, такие как, методы скользящей средней, экспоненциального сглаживания и др. ^[9], в частности методы, связанные с добавлением белого шума и последующим устранением шумовой компоненты за счет его усреднения ^[10].

2. Оценка шумовой компоненты в исходных данных

Обозначим значения аргумента (временные отрезки) через xi, соответствующие им наблюдаемые значения через y_i, предполагаемую функцию «истинной» зависимости – f(x). И исследуемый ряд представим в виде пар (x_i , f(x_i)). Будем рассматривать общий случай, когда исходные данные задаются с переменным шагом x_i+1=x_i+h_i, где величины h_i, вообще говоря, различны.

Наблюдаемые данные будем представлять в виде , где s_i – шум.

Пусть функцию f(x) имеет производные до 4 порядка включительно, тогда ее значение в точке x+t может быть представлено как

(1)

Рассмотрим значения функции в окрестности точки x₀. Выберем точки x_k=x₀+t_k, где k=0,1,2,3,4; t₀=0 и все t_k различны.

В матричном виде она примет вид.

AB=C

Где ; ;.

Индекс k принимает значения 1, 2, 3, 4.

Ранг матрицы A равен 3, следовательно, существует вектор , где не все λ_k равны 0 такой, что , тогда .

Значения λ_k определяются с точностью до постоянного множителя, в частности, положив λ₄=1, получим допустимый набор из решения системы линейных уравнений

(2)

Определитель матрицы A

Аналогично

Поскольку вектор определен с точностью до произвольного множителя, то их удобнее представить в симметричном виде

(3)

В частности, решение этой системы для равноотстоящих узлов ^[11] t₀=0, t₁=-1, t₂=1, t₃=-2, t₄=2 дает значения

Учитывая (2) получаем

(4)

Отметим, что выражение слева дает с точностью до множителя t₄ численное представление четвертой производной и при функциях f(x), представимых в виде полиномов не выше третьей степени, тождественно обращается в 0.

Далее учтем, что f(x_m+t_k)=y_m+k-s_m+k и, следовательно

При отсутствии быстрых, то есть с периодами соизмеримыми с шагом ряда, осцилляций f(x) величину можно считать малой, и тогда получаем

(5)

Величина систематических отклонений от 0 будет тем меньше, чем меньше будет , то есть при выборе набора из пяти точек в качестве базовой точки x целесообразно выбирать точку x₃. Перенумеруем точки и введем , тогда (5) примет вид

(6)

Правая часть равенства (6) представляет собой случайную величину. Пусть величины s_m+k – независимые случайные величины с 0 математическим ожиданием и дисперсией σ², тогда математическое ожидание .

Среднеквадратичное значение шума σ2, таким образом, можно оценить исходя из

(7)

3. Выделение в данных функциональной и шумовой компонент

В этих условиях можно определить значения шумовой компоненты s и, следовательно, f(x) исходя из

(8)

Перепишем (7) в матричных обозначениях. Для этого введем матрицу L=(l_ij), где , тогда

(9)

И условие минимума принимает вид

или

(10)

В отличие от случая с равноотстоящими узлами здесь возникает необходимость вычисления элементов матрицы L для каждой строки (при равноотстоящих узлах ненулевые элементы матрицы L одинаковы для всех строк). Система (9) из-за ограничения является нелинейной и нахождение ее решений даже при небольшой размерности встречает значительные вычислительные трудности. В то же время решение системы линейных алгебраических уравнений (СЛАУ) из (10) при заданном значении t не создает особых проблем. В самом деле, матрица является 9 диагональной ленточной матрицей и можно показать, что время решения СЛАУ с такой матрицей линейно по числу уравнений.

Отметим, что матрица является неотрицательно определенной и все ее собственные числа μ_i≥0.

Пусть E_i – собственные вектора матрицы, соответствующие собственным числам μ_i≥0,

Тогда первое из равенств (9) примет вид

В этих условиях представляет собой при t>0 монотонно убывающую функцию от t.

Заметим, что при t→0 минимум функции F(s) в (8) будет достигаться, как следует из (5), при . Последнее условие означает стремление 4-ой производной к 0, то есть исходная функция будет близка к полиному не выше 3 степени.

При t→∞ минимум F(s) будет достигаться, при s_i→0 и, следовательно, функция f(x_i)→y_i.

4. Алгоритм удаления шума и выделения функциональной компоненты в данных

Учитывая сделанные замечания о характере функции F(s) можно предложить следующий итеративный алгоритм.

1. Вычисляем оценку шума σ² на основе (6). Задаем начальное значение t₀.

2. Решаем СЛАУ в соответствии с (9) и вычисляем значение <s,s>.

3. Если <s,s> > σ², переходим к пункту 4, иначе к пункту 5.

4. Решаем СЛАУ в соответствии с (9) с t=2t₀ и вычисляем значение <s,s>. Если <s,s> >nσ², устанавливаем t₀=t и повторяем пункт 4. В противном случае устанавливаем t₁=t и переходим к пункту 6.

5. Решаем СЛАУ в соответствии с (9) с t=t₀/2 и вычисляем значение <s,s>. Если <s,s> <nσ², устанавливаем t₀=t и повторяем пункт 5. В противном случае устанавливаем t₀=t, t₁=t₀ и переходим к пункту 6.

6. Искомое значение t лежит между t₀ и t₁. Решаем СЛАУ в соответствии с (9) с t= (t₀ + t₁)/2 и вычисляем значение <s,s>. Если |<s,s> -nσ²|<ε, то останавливаем процесс и на основе вычисленных значений s_i рассчитываем значения искомой функции f(x_i)= y_i-s_i. Если нет, то строим новый интервал, устанавливая в зависимости от выполнения неравенства <s,s> >nσ² либо t₀=t, либоt₁=t и переходим к пункту 6.

Отметим, что говоря о шуме и его дисперсии в (6) и (7) предполагается, что он представляет собой одинаково распределенную случайную величину на всем временном ряде. Если это не так, то вместо абсолютной величины шума s_i в (5) шум удобнее представлять в виде s_i=φ(x,y)u_i так, чтобы случайная величина u_i была бы одинаково распределенной на всем временном ряде. Наиболее естественно предполагать при большом разбросе исходных данных, что шум является результатом измерений, которые имеют постоянной относительную погрешность. В этом случае естественно принять s_i=y_iu_i, тогда (6) и (7) примут соответственно вид

(10)

(11)

Тогда дисперсия относительного шума u определяется из (11)

И задача (7) приобретает вид

(12)

Или в матричном виде

Тогда условие минимума принимает вид

или

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

(13)

В частности, решение этой системы для равноотстоящих узлов ^[11] t₀=0, t₁=-h, t₂=h, t₃=-2h, t₄=2h дает значения

5. Результаты численного моделирования

На рисунке 1 представлены результаты обработки зашумленных данных.

Значения по осям x и y – , где rnd представляет собой случайную величину, равномерно распределенную на интервале (-0,1;0,1).

Рисунок 1

На рисунке 2 представлены результаты обработки данных, представляющих значения функции y=e^x на интервале (0;10) с шагом 0,1 округленные до двух значащих цифр.

Среднеквадратичная абсолютная погрешность исходных данных от теоретической кривой в рассматриваемом примере составляет 56, при удалении абсолютного шума - 30, при удалении относительного шума – 19. Таким образом, в случаях, когда данные меняются в широких пределах, а в данном случае отношение максимального значения к минимальному составляет 22000, удаление относительного шума приводит и к лучшему удалению абсолютного шума по сравнению с методом прямого удаления абсолютного шума.

Рисунок 2

6. Выводы

Таким образом, предлагаемая методика и алгоритмы выявления и устранения шума в данных в предположении о гладкости, представляемой ими функции, позволяют:

· обоснованно определить уровень как абсолютного, так и относительного шума в данных вне зависимости от равномерности шага измерений в исходных данных;

· удалить из данных шумовую компоненту;

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