6.15M

Программный контроль временных рядов на наличие выбросов

1.

Лабораторная работа № 2
Программный контроль
временных рядов
на наличие выбросов
Профессор Кузнецов Анатолий Дмитриевич
Доцент Сероухова Ольга Станиславовна
Доцент Восканян Карина Левановна
Учебная дисциплина
«АВТОМАТИЧЕСКИЕ МЕТЕОРОЛОГИЧЕСКИЕ СТАНЦИИ ОБЩЕГО И СПЕЦИАЛЬНОГО НАЗНАЧЕНИЯ»
2019

2. 3. Цель работы:

1.
Произвести программный контроль на наличие
выбросов во временных рядах метеорологических
величин, измеренных с помощью АМС.
2.
Проанализировать полученные результаты.

4.

Порядок выполнения:
1. Сформировать таблицу с требующими проверки
временными рядами в кодировке «Exсel» (в случае
необходимости переведя их из текстового формата в
таблицу «Exсel»).
Рис. 1. Пример представления временных рядов, сформированных
для последующего контроля на выбросы

5.

2. Запустить файл «Контроль 4М по отрезкам (Ирвин
1).xlc» и стереть всю имеющуюся здесь информацию на
Листах 1 и 2 (до выполнения этой операции необходимо
ознакомиться с формой представления данных
расчетов на Листах 1 и 2).
В том числе удалить и имеющиеся на Листе 2 все
графиками – их в дальнейшем надо будет заново
построить с учетом собственных проверяемых данных!
Рис. 2. Пример фрагмента Листа 2 после запуска файла «Контроль 4М по
отрезкам (Ирвин 1).xlc»

6.

3. Переписать требующие проверки временные ряды из
ранее созданного файла (см. рис. 1) на Лист 1 файла
«Контроль 4М по отрезкам (Ирвин 1).xlc». В
дальнейшем программой может быть осуществлена
проверка каждого из этих рядов по отдельности (каждая
колонка данных отдельно).
4. Определить на Листе 1 номер колонки с подлежащим
проверки временным рядом (колонке «А» соответствует
номер 1, колонке «В» соответствует номер 2 и т.д., )
5. Определить общее число элементов во временном
ряде, подлежащим проверке. Например, N = 7712 .

7.

6. Запустить макросы.
В диалоговом режиме ввести все запрашиваемые
программой
параметры
(при
каждом
запуске
указывается только один номер колонки).
Всего необходимо ввести 6 параметра:
- номер колонки с проверяемым рядом на Листе 1;
- общую длину проверяемого ряда N;
- начальное значение, с которого начнется проверка
ряда: n1;
- конечное значение,которым закончится проверка
ряда: n2;
- длину отрезка временного ряда, для которого
последовательно будут вычисляться среднее значение
и среднеквадратическое отклонение (СКО);
- коэффициент кско, позволяющий рассчитать
критическое значение СКО отрезка временного ряда:
если фактическое значение СКО отрезка временного
ряда в кско раз меньше СКО для всего ряда, то оно не
используется для контроля за выбросом.

8.

Этап 1.
Контроль
на наличие пропусков

9.

После
ввода
последнего
параметра
программой осуществляется контроль наличия
пропусков в анализируемом ряду.
При обнаружении хотя бы одного пропуска
дальнейшая работа по контролю ряда
прекращается и может быть проведена только
после устранения наличия пробелов во
временном ряду (см. рис. 3 и 4).

10.

Рис. 3. Пример вида Листа 1 после работы программы в том случае, когда
во временном ряде присутствуют разрывы
(в данном случае под номером 644 – см. колонку «I»)

11.

Рис. 4. Пример вида Листа 2 после работы программы в том случае, когда во
временном ряде присутствуют разрывы
(до устранения разрывов дальнейший контроль не производится)

12.

Если пробел соответствует небольшому временному
промежутку, то он может быть заполнен либо на основе
временной
интерполяции,
либо
другими
существующими способами (см. рис. 5 и 6).

13.

а)
б)
Рис. 5. Пример заполнения пропуска во временных рядах
(а – фрагмент исходных рядов, б – после заполнения пропуска путем
удаления пустой строки)

14.

а)
б)
Рис. 6. Пример заполнения пропуска во временных рядах
(а – фрагмент исходных рядов, б – после заполнения пропуска путем
линейной интерполяции)

15.

Этап 2.
Контроль
на наличие «выбросов»

16.

В программе предусмотрено четыре метода
обнаружения
резких
изменений
(«выбросов»)
метеорологического параметра во временном ряде.
Каждый из них имеет свои особенности, что
приводит к весьма большому разнообразию в
полученных результатах. Последнее делает анализ
полученных результатов весьма не простой задачей.
Примером разнообразив полученных результатов
может служить рис. 7, где, в частности, на 4 графиках
представлены номера подозрительных на наличие
выбросов отсчетов во временном ряде температуры,
полученные при использовании четырех методов
контроля ряда на наличие «выбросов».

17.

Рис. 7. Пример вида Листа 2 после работы программы

18.

Окончательное решение о том, является ли то или
иное определенное программой значение временного
ряда «выбросом» остается за исследователем!
В частности, в этом помогает комплексный контроль
синхронных
временных
рядов
различных
метеорологических величин (например, температуры и
давления), построение графиков отрезков временных
рядов в районе предполагаемого выброса и др.

19.

Метод 1
Сравнение
двух значений ряда
на границе
контрольного участка
(с использование СКО
контрольного участка)

20.

При рассмотрении этого метода учтем, что при
работе программы для анализа выбросов может
быть выбран фрагмент из всего ряда, ограниченного
значениями n1 и n2. В дальнейшем, при
рассмотрении других методов, будем считать, что
анализируется весь временной ряд.
При использовании этого метода на участке
временного ряда от n1 до n2 используются
введенные ранее:
- длина «контрольного участка» временного ряда
ndl ;
- коэффициент кско.

21.

Алгоритм метода 1:
1. На участке от n = n1 до n
рассчитывается величина СКО.
=
n1+ndl-1
2. Затем последнее значение этого контрольного
участка (i = n1+ndl-1 ) сравнивается с граничным для
этого отрезка значением (i = n1+ndl ) и по абсолютной
величине их разности определяется, во сколько раз она
превосходит среднее квадратическое отклонение (СКО):
от 3 до 6 раз.
При этом, если вычисленное значение СКО (от n = n1
до n = n1+ndl-1 ) оказывается в кско раз меньше СКО
для всего ряда (от n = n1 до n = n2), то сравнение
абсолютных разностей не производится - на
контрольном участке изменчивость слишком мала и не
может служить критерием наличия или отсутствия
выброса.

22.

4. Затем границы контрольного участка сдвигаются
на один шаг вперед и процедура оценивания
повторяется: рассматривается отрезок ряда от n = n1+1
до n = n1+ndl .
5. Процесс проверки завершается, когда граничным
значением для участка длинной ndl становится
последний элемент проверяемого фрагмента ряда:
n = n2.

23.

12.6
Следующий контрольный
участок протяженностью ndl
12.4
12.2
асл
12
всл
11.8
11.6
11.4
11.2
a
b
11
10.8
10.6
0
5
10
15
20
25
30
35
Контрольный участок
протяженностью в ndl
отсчетов: расчет СКОndl
Рис. 8. Иллюстрация схемы расчета критерия наличия выброса по методу 1:
критерий Р = [Abs(a - b)] / (СКОndl)
(а – последний отсчет в контрольном участке, в – проверяемое на наличие
выброса значение временного ряда; асл и всл – отсчеты, используемые на
следующем шаге контроля)

24.

Анализ данных с использованием метода 1.
1. Для визуального контроля проверяемого временного
ряда строим его график (см. рис. 9). На рис. 9 мы видим
несколько участков, «подозрительных» на наличие
выбросов. На рисунке они отмечены стрелками.
15
Т, С
10
5
0
-5
1
1000
1999
2998
3997
4996
5995
Порядковый номер измерения
Рис. 9. График проверяемого временного ряда
6994

25.

Для предварительного определения порядкового номера того
или иного «подозрительного» участка наведите на соответствующую
точку курсор – высветится порядковый номер этой точки (см. рис. 10)
Рис. 10. Определение порядкового номера «подозрительной» точки

26.

Индикатор выброса (СКО)
2. Следующий этап анализа – построение графика данных,
содержащихся в колонке 2 на Листе 2 (рис. 11). По горизонтальной
оси – порядковый номер измерения, по вертикальной оси –
индикатор выброса Р. Численное значение Р показывает, во сколько
раз абсолютная разность последнего значения контрольного участка
и следующего значения превосходит величину СКО, вычисленного
для этого контрольного участка. В колонке «D» и «Е» на Листе 2
можно посмотреть порядковые номера измерений, для которых
величина Р превосходит значения 6 и 5 соответственно (см. рис. 12).
7
6
5
4
3
2
1
0
1
1000
1999
2998
3997
4996
5995
6994
Порядковый номер измерения
Рис. 11. График индикатора наличия
«подозрительных» точек
Рис. 12. Данные контроля методом 1.

27.

3. С учетом полученных данных анализируем
поведение ряда в области выявленной точки с номером
3729: анализируем имеющиеся в этой области значения
ряда (см. рис. 13) и строим вспомогательный график
имеющейся в этой области временной изменчивости ряда
(см. рис. 14).
По этим данным, с учетом физического анализа – т.е.
с учетом того, что исследуемый временной ряд
представляет собой значения приземной температуры с
дискретностью 5 мин, делаем вывод, что в ряде
присутствует
«псевдовыброс»,
связанный,
по
всей
видимости, с временным разрывом в работе АМС.

28.

Следовательно, по этим данным данный ряд можно
считать однородным:
с n = 1 по n = 3728
и
с n = 3729 по n = 7712.
3
2.5
2
1.5
1
0.5
0
-0.5
-1
-1.5
Рис. 13. Фрагмент
временного ряда в
районе
«подозрительной» точки
3725
3726
3727
3728
3729
3730
3731
3732
3733
Рис. 14. График фрагмента временного
ряда в районе «подозрительной» точки

29.

Метод может давать ложную «тревогу» на участках, где
наблюдается практически постоянные значения метеорологической
величины. Проиллюстрируем это на следующем примере.
После работы программы был получен следующий график
данных, содержащихся в колонке «В»: номер подозрительного на
выброс измерения равен 2254, разность двух соседних измерений в 3
раза превосходит СКО контрольного участка.

30.

Рассмотрим фрагмент временного ряда в районе
подозрительной на выброс точки: 2254.
-6.04
-6.06
-6.08
Т, С
-6.1
-6.12
-6.14
-6.16
-6.18
-6.2
-6.22
2245
2246
2247
2248
2249
2250
2251
2252
2253
2254
2255
2256
2257
2258
2259
2260
Порядковый номер измерения
График фрагмента временного ряда в районе
«подозрительной» точки (диапазон изменения
температурной шкалы: от -6.04 до -6.22, Δtшкалы = 0.18 0С)
-5
-5.2
-5.4
-5.6
-5.8
Т, С
Фрагмент временного
ряда в районе
«подозрительной» точки
2254 на Листе 2,
изменение температуры
за 15 мин равно 0.1 0С
-6
-6.2
-6.4
2245
2246
2247
2248
2249
2250
2251
2252
2253
2254
2255
2256
2257
2258
2259
2260
-6.6
-6.8
-7
Порядковый номер измерения
График фактического изменения температуры в
районе «подозрительной» точки (диапазон
изменения т-ры: от -6.0 до -7.0, , Δtшкалы = 1 0С)

31.

Метод 2
Сравнение
последнего в контрольном
участке значения ряда
со средним значением
для этого контрольного участка
(с использование СКО
контрольного участка)

32.

При реализации этого метода используются
введенные ранее:
- длина «контрольного участка» временного ряда
ndl ;
- коэффициент кско.

33.

Алгоритм метода 2:
1. Для участка временного ряда с номерами от n = 1
до n = ndl рассчитываются: среднее значения Srndl и
величина СКОndl.
2. Затем граничное значение этого контрольного
участка (n=ndl+1) сравнивается со средним для этого
отрезка значением Srndl и по абсолютной величине их
разности определяется, во сколько раз она превосходит
среднее квадратическое отклонение (СКОndl): от 3 до 6
раз.
При этом, если вычисленное значение СКОndl
оказывается в кско раз меньше СКОN для всего ряда (от
n = 1 до n = N) , то сравнение абсолютных разностей не
производится - на контрольном участке изменчивость
слишком мала и не может служить критерием наличия
или отсутствия выброса.

34.

3. Затем границы контрольного участка сдвигаются
на один шаг вперед и процедура оценивания
повторяется: рассматривается отрезок ряда от
n = 2 до n = ndl +1.
4. Процесс проверки завершается, когда граничным
значением для участка длинной ndl становится
последний элемент ряда:
n = N.

35.

12.6
Следующий контрольный
участок протяженностью ndl
12.4
асл
12.2
12
11.8
11.6
11.4
11.2
a
11
10.8
10.6
0
5
10
15
20
25
30
35
Контрольный участок
протяженностью в ndl
отсчетов: расчет среднего
Srndl и СКОndl
Рис. 15. Иллюстрация схемы расчета критерия наличия выброса по методу 2:
критерий Р = [Abs(Srndl - a)] / (СКОndl)
( а – проверяемое на наличие выброса значение временного ряда; асл – отсчет,
используемый на следующем шаге контроля)

36.

Анализ данных с использованием метода 2.
1. Для визуального контроля проверяемого временного
ряда строим его график (см. рис. 16). На рис. 15 мы
видим несколько участков, «подозрительных» на наличие
выбросов. На рисунке они отмечены стрелками.
15
Т, С
10
5
0
-5
1
1000
1999
2998
3997
4996
5995
Порядковый номер измерения
Рис. 16. График проверяемого временного ряда
6994

37.

2. Следующий этап анализа – построение графика данных,
содержащихся в колонке 3 на Листе 2 (рис. 17).
По горизонтальной оси – порядковый номер измерения, по
вертикальной оси – индикатор выброса Р.
Численное значение Р показывает, во сколько раз
абсолютная разность последнего значения контрольного участка и
следующего значения превосходит величину СКО, вычисленного для
этого контрольного участка.
В колонке «F» и «G» на Листе 2 можно посмотреть
порядковые номера измерений, для которых величина Р
превосходит значения 6 и 5 соответственно (см. рис. 18).

38.

7
Индикатор выбросов (ско)
6
5
4
3
2
1
0
1
1000
1999
2998
3997
4996
5995
6994
Порядковый номер измерения
Рис. 17. График индикатора наличия
«подозрительных» точек
Рис. 18. Данные контроля методом 2.

39.

Отмечаем, что в отличие от метода 1, в данном случае
выявлено существенное большее количество точек, подозрительных
на «выброс». Для Р = 6: 1 точка, для Р = 5: 3 точки, для Р = 4: 10
точек и для Р = 3: 44 точки.
При этом следует учитывать, что при подсчете количества
точек, например, для Р = 4, учитывались точки с Р = 4 плюс
количество точек с Р = 5 и плюс количество точек с Р = 6.
Таким образом можно сделать вывод о том, что метод 2
более чувствителен к изменчивости временного ряда, поскольку
выявил значительно большее число «подозрительных» точек.
Однако напомним, что как и для метода 1, окончательное решение о
том, является ли то или иное определенное программой значение
временного ряда выбросом остается за исследователем!
В дальнейшем остановимся на исследовании поведения
ряда в районе «подозрительных» точек, для которых Р = 5 и Р = 6.

40.

3. С учетом полученных данных анализируем поведение ряда в
области выявленной точки с номером 3729 (значение индикатора Р = 6),
и двух точек с номерами 3239 и 3730 (значение индикатора Р = 5). Для
этого анализируем имеющиеся в этой области значения ряда (см. рис.
19) и строим вспомогательные графики имеющейся в этих областях
временной изменчивости ряда (см. рис. 20).
3.5
3
2.5
а)
2
1.5
1
0.5
0
3234
3233
3232
3231
3230
3229
3228
3227
3226
3225
3
а)
б)
Рис. 19. Фрагменты временного
ряда в районе «подозрительных»
точек
2.5
2
1.5
б)
1
0.5
0
-0.5
-1
-1.5
3725
3726
3727
3728
3729
3730
3731
3732
3733
Рис. 20. Графики фрагментов временного
ряда в районе «подозрительных» точек

41.

По этим данным, с учетом физического анализа – т.е. с учетом
того, что исследуемый временной ряд представляет собой значения
приземной температуры с дискретностью 5 мин, делаем вывод, что в
ряде между отсчетами n = 3728 и n = 3729 присутствует «выброс»,
связанный, по всей видимости, с временным разрывом в работе АМС
(см. рис. 17 б и рис. 18 б). Этот вывод совпадает с тем, который был
получен при использовании метода 1.
Природа
появления
при
использовании
метода
2
«подозрительной» точки с номером 3730 (Р = 5) так же связана с этим
фактом – это «отголосок» наличия большого градиента между точками
3728 и 3729 .
Для анализа «подозрительной» точки с номером 3239 обратимся
к рис. 19а и 29а. В этом случае вывод о наличии сбоя в работе АМС в
районе номера 3239 уже не так однозначен и для окончательного
решения о характере этой точки необходим дополнительный анализ.

42.

Метод 3
Метод 1-3 сигм
для контрольного участка
(с использование среднего и СКО
контрольного участка)

43.

Вероятности того, что отклонение распределенной по
нормальному закону непрерывной случайной величины X, от
математического ожидания a по абсолютной величине (то есть по
модулю) будет меньше δ , определяется следующей формулой:
где a - математическое ожидание, σ - среднее квадратическое
отклонение.
Табличные значения функции Ф(x):
Ф(1) = 0.34134; Ф(2) = 0.47725; Ф(3) = 0.49865/
Следовательно, вероятность того, что случайная величина,
имеющая нормальный закон распределения, отклонится от своего
математического ожидания на величину, определяемую кратностью
среднего квадратичного отклонения, равна:
±1·СКО => Ф(1)·2 = 68.3%;
±2·СКО => Ф(2)·2 = 95.5%
±3·СКО => Ф(3)·2 = 99.7%
Откуда следует, что вероятность того, что случайная величина
отклонится от своего математического ожидания на большую
величину, чем утроенное среднее квадратичное отклонение,
практически равна нулю.

44.

К сожалению, временные ряды метеорологических величин не
имеют нормального распределения, хотя в ряде случаев близко к
нему.
Тем не менее, для оценки выбросов можно воспользоваться
оценками, основанными на использовании абсолютных отклонений
значений временного ряда от среднего значения.
Предполагая, что с вероятностями, близкими к значениям 68%;
95% и 99% те значения временного ряда, которые отличаются от
среднего значения Srndl в 1, 2 или 3 раза от СКОndl, являются
выбросами – т.е. не типичными значениями для рассматриваемой
выборки протяженностью ndl отсчетов.
На таком предположении и основан третий из рассматриваемых
методов выявления выбросов в исследуемом временном ряде.

45.

12.6
a2
12.4
a3
12.2
andl
12
Следующий контрольный
участок протяженностью ndl
a1
11.8
11.6
11.4
11.2
11
10.8
10.6
0
5
10
15
20
25
30
35
Контрольный участок
протяженностью в ndl
отсчетов: расчет среднего
Srndl и СКОndl
Рис. 21. Иллюстрация схемы расчета критерия наличия выброса по методу 3:
критерий Р = [Abs(Srndl - ai)] / (СКОndl), i = 1, 2, … , ndl
( аi – проверяемое на наличие выброса значение временного ряда)

46.

Метод Ирвина
Сравнение
крайних значений
вариационного ряда
для контрольного участка
с критерием Ирвина
(с использование СКО
контрольного участка)

47.

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

48.

На рис. 22 представлен исходный фрагмент (выборка)
временного ряда (черная линия), содержащий 50 значения
температуры воздуха (отсчет под номером 43 является
«выбросом») и соответствующий ему вариационный ряд (красная
линия), построенный в порядке возрастания.
20
19
18
T, C
17
Ряд температуры
16
Вариационный ряд
15
14
13
12
1
3
5
7
9 11 13 15 17 19 21 23 25 27 29 31 33 35 37 39 41 43 45 47 49
Порядковый номер
Рис. 22. Исходный и вариационный ряды температуры воздуха
(пояснения в тексте)

49.

В вариационном ряду значений выборки оценивают
сомнительное значение на одном из концов ряда. Для этого
вычисляют расчётное значение критерия Ирвина Iрасч для
максимального значения вариационного ряда:
Iрасч= Abs(хndl – хndl – 1 )/СКОndl,
и для минимального значения вариационного ряда:
Iрасч= Abs(х1 – х2 )/СКОndl,
где х1 и хndl – сомнительные значения,
х2 и хndl - 1 –значения в вариационного ряда, соседствующие с
сомнительными значениями;
СКОndl - выборочное среднеквадратическое отклонение,
рассчитываемое для контрольного участка протяженностью ndl с
учётом сомнительного значения.

50.

Полученное расчётное значение сравнивают с табличным Iтабл ,
которое зависит от объёма выборки n и принятого уровня значимости α.
Если Iрасч > Iтабл, то сомнительное значение считают грубой
ошибкой.
Однако при автоматизированной обработке данных удобно
рассчитывать Iтабл с приемлемой точностью по соотношениям,
показанным в табл. для случая использования выборочного СКО
при изменении объёма выборки n в пределах от 3 до 1000.
α
0,01
0,05
0,1
Iтабл
-205,06n-3 + 424,26n-2,5 - 352,483n-2 +143,747n-1,5 - 33,401n-1+6,381n-0,5 + 1,049
-229,21n-3 + 422,39n-2,5 - 320,96n-2 +124,594n-1,5 - 26,15n-1+4,799n-0,5 + 0,7029
-132,78n-3 + 224,24n-2,5 - 165,27n-2 +68,614n-1,5 - 16,109n-1+3,693n-0,5 + 0,549