Идентификация тренда временного ряда

Лекция 11 Идентификация тренда временного ряда Прежде, чем искать причину какого-либо факта, необходимо как следует убедиться в том, что такой факт действительно имеет место. Фонтенель «История оракулов»,1687 Идентификация тренда ВР связана с нахождением аналитического вида его математической модели. При этом принимается допущение, что «ошибки» измерений, в силу рассмотренных ранее причин, носят случайный характер и распределены примерно по нормальному закону с нулевым математическим ожиданием и с некоторой (заранее неизвестной) дисперсией  , а сами измерения подчиняются равномерному (равновероятному) закону в некоторых заранее известных пределах (рис.36). 2 X X max РЗ X cp НЗ X min t t1 tm Рис.36 Поставим задачу идентификации тренда ВР в виде аппроксимирующей функции X~ (t ) , которая «в среднем совпадает» с результатами измерений в области значений Xj и tj,, где j=1,…,m. Заметим, что в такой постановке задача аппроксимации отличается от задачи интерполяции. При интерполяции известны точные значения функции в n узловых точках. Поэтому интерполирующая функция, например, в форме полинома степени (n-1), должна точно совпадать с заданными точными значениями во всех n узловых точках. Аппроксимирующая функция, в отличие от интерполирующей функции, не проходит точно через узловые точки ВР, она проходит «в среднем» через поле точек и представляет собой характеристику «истинного» сигнала, «отфильтрованного» от шумов.Последнее подробно будет рассмотрено несколько позже. Как было определено, ВР есть множество {𝑋𝑗 }𝑚 𝑗=1 из m чисел, представляющих значения выборок некоторого измеренного параметра 𝑋(𝑡𝑗 ). Такое множество, будучи представленным на плоскости (X,t), как правило, характеризуется некоторым разбросом точек с соответствующими координатами, по которым трудно выдвинуть гипотезу о подходящем типе функции X(t), что можно проиллюстрировать рис.37,а. Назовём множество { X j , t j } ансамблем A(X,t), если такое предположение осуществляется. Ансамбль (фр. ensemble — совокупность, согласованность, единство частей, образующих что-либо целое). В задаче прогнозирования технического состояния КЭО установить наличие ансамбля в ВР можно, если удаётся найти такую фазовую плоскость (𝜂, 𝜇), где 𝜂 = 𝜂(𝑋, 𝑡) и 𝜇 = 𝜇(𝑋, 𝑡), на которой точки с координатами 𝜂𝑗 и 𝜇𝑗 с некоторым допустимым разбросом располагаются вдоль прямой или некоторой гладкой кривой (рис.37,б).   ( X ,t) а б С   t ( X , t) Рис.37 Когда установлено наличие ансамбля, можно приступать к решению задачи идентификации функции 𝜂(𝜇). Затем, действуя методом «обратного хода», можно оценить тип и параметры искомой функции тренда 𝑋̃(𝑡) – аппроксимирующую функцию BP. При этом важно, чтобы функции 𝜂(𝑋, 𝑡) и 𝜇(𝑋, 𝑡) обеспечивали не только принципиальную возможность перехода от 𝜂(𝜇) к X(t), но и осуществляли это, желательно, наиболее простыми способами. Далее будет рассмотрен уже знакомый вам материал применительно к задаче прогнозирования. Поскольку тренд ВР, подлежащий идентификации, отличается разбросом точек на плоскости (Х,t), то естественно применить какую-то операцию сглаживания, например, интегрирование Ik, где k – кратность интегрирования. Интеграл кратности k берётся в пределах от 0 до t (или в других пределах). Следует особо отметить, что операция интегрирования, по сути, есть преобразование 𝐼𝑘 (𝑋, 𝑡) → 𝜇(𝑋, 𝑡), а исходный ряд есть 𝜂(𝑋, 𝑡). Рассмотрим пример. На рис.38 показан результат процедуры сглаживания (Хс) некоторого временного ряда X(t) оператором 𝐼 c k=1. X,Xc Хс 140 6 Х 120 5 100 4 80 3 60 2 40 1 20 1 2 3 4 5 6 7 8 t 9 10 11 12 13 14 15 Рис.38 Как видно, исходный ВР (чёрная ломаная линия Х, образованная с учётом разброса точек на плоскости) отличается существенно неравнозначными по значениям выборками. Сглаженный ВР (красная пунктирная линия) уже не имеет существенных различий в значениях, он практически монотонен. На рис.39 показана фазовая плоскость 𝜂 = (Х, 𝑡), 𝜇 = 𝐼(Х, 𝑡) с расположением точек с координатами Х𝑗 , 𝐼(𝑋𝑗 ), 𝑗 = 1,2, . . . , 𝑚. На рис.39 видно, что поле точек образует ансамбль, т.к. точки группируются с некоторым рассеянием вдоль прямой CD, которую можно представить аналитически так: 𝑡 𝑋(𝑡) = а0 + а1 𝐼[𝑋(𝑡)] = а0 + а1 Х̄ (𝑡) = а0 + а1 ∫ 𝑋(𝜏)𝑑𝜏. (3.1) Заметим, что здесь не указывается кратность интеграла.  D 6 5 4 3 2 1 С 20 40 60 80 100 120 140 160  Рис.39 В выражении (3.1) 𝑋̄(𝑡) - более краткое обозначение операции Iинтегрирования; а0 – параметр, характеризующий начальные условия процесса деградации технического состояния (работоспособности) КЭО или его узла, например, начальный уровень вибрации подшипникового узла или косвенный параметр, определяющий зазор между элементами узла КЭО и т.п. Параметр а1 (наклон прямой) определяет темп развития деградации на временном интервале индикации диагностического параметра (в рассматриваемом примере 𝑡 = [0,15[ ) и называется параметром тренда. Выражение (3.1) есть интегральное уравнение, общее решение которого имеет вид: Х(𝑡) = 𝑎0 𝑒𝑥𝑝( 𝑎1 𝑡) = Х̃ (𝑡) . (3.2) Выражение (3.2) и есть результат «обратного хода», позволяющий определиться с типом функции тренда ВР. Эта функция тренда является аппроксимирующей функцией 𝑋̃(𝑡) исходного ВР контролируемого (индицируемого) параметра. Свойствами сглаживания кроме обычного классического I-интеграла обладает и L –интеграл (см. ниже справку): 𝑡 𝐿[𝑋(𝑡)] = Х̑ (𝑡) = ∫0 Х(𝜏)𝑑𝜏 + 𝑡 ⋅ Х(𝑡) (3.3) В L - интеграле первый суммант (обычный интеграл) отражает как бы всю сглаженную историю тренда ВР Х(t), а второй суммант tX(t) характеризует значимость (вес) каждой выборки ВР. Как видно, значимость конкретного члена ВР определяется не только числовым значением, но положением на оси времени, что справедливо для процесса деградации технического состояния объекта, т.к. последние значения ВР более объективно характеризуют состояние КЭО (его контролируемого узла). Справка. Использовать L-интеграл как оператор сглаживания в задачах нахождения экспериментальных трендов впервые предложил заведующий кафедрой Судовой автоматики и измерений ЛКИ, профессор Соболев Л.Г. (1926 - 2007), разработавший на базе обобщенного операционного исчисления теорию L-преобразований. Более подробно с теорией L-преобразований можно ознакомиться, например, в учебном пособии “Обработка результатов измерений в судостроении”, изд. ЛКИ, 1983 г. Линейная зависимость в фазовой плоскости представляется в развернутом виде уравнением: 𝑋(𝑡) = 𝑎0 + 𝑎1 𝐿[𝑋(𝑡)] = 𝑎0 + 𝑎1 ⋅ 𝑋̑(𝑡) = 𝜂(𝑋, 𝑡); 𝜇[𝐿(𝑋, 𝑡)] 𝑡 = 𝑎0 + 𝑎1 ⋅ [∫0 𝑋(𝜏)𝑑𝜏 + 𝑡 ⋅ 𝑋(𝑡)]. Решение этого уравнения в общем виде есть выражение 𝑋(𝑡) = 𝑎0 /(1 − 𝑎1 ⋅ 𝑡)2 = 𝑋̃(𝑡) (3.4) (3.5) и также является альтернативной аппроксимирующей функцией Х̃ (𝑡) наблюдаемого ВР диагностического параметра по отношению к выше рассмотренной. При 𝑡 = 1⁄𝑎1 тренд- функция (3.5) принимает значение 𝑋̃(t)→ ∞, что можно трактовать как параметрический отказ КЭО. Численные значения параметров аппроксимирующих определяются известным методом наименьших квадратов: функций 2 𝑚 2 𝑚 𝑚 𝑚 𝑎0 = [∑𝑚 𝑗=1 𝑋с 𝑗 ∑𝑗=1 𝑋𝑗 − ∑𝑗=1 𝑋с 𝑗 ∑𝑗=1 𝑋с 𝑗 𝑋𝑗 ]/[𝑚 ∑𝑗=1 𝑋с 𝑗 − 2 (∑𝑚 𝑗=1 𝑋с 𝑗 ) ]; 2 𝑚 𝑚 𝑚 𝑚 2 𝑎1 = [𝑚 ∑𝑚 𝑗=1 𝑋с 𝑗 𝑋𝑗 − ∑𝑗=1 𝑋с 𝑗 ∑𝑗=1 𝑋𝑗 ]/[𝑚 ∑𝑗=1 𝑋с 𝑗 − (∑𝑗=1 𝑋с 𝑗 ) ], (3.6) где Х𝑐𝑗 = 𝐼[𝑋(𝑡𝑗 )] ∨ 𝐿[𝑋(𝑡𝑗 )].