CC BY
107
Изв. вузов «ПНД», т. 17, № 5, 2009
УДК 530.182+001.891.57
МОДЕЛИРОВАНИЕ ПО ВРЕМЕННЫМ РЯДАМ В ПРИЛОЖЕНИИ К ОБРАБОТКЕ СЛОЖНЫХ СИГНАЛОВ
Б.П. Безручко
Сигналы, снимаемые с большинства реальных систем, особенно живых, нерегулярны, часто хаотичны, нестационарны и зашумлены. Так как современные измерительные приборы обычно реализуют цифровую обработку информации, записи сигналов имеют вид дискретной последовательности отсчетов (временных рядов). В данной лекции дается краткий обзор возможностей обработки таких экспериментальных данных с опорой на реконструкцию и использование предсказательной модели рассматриваемой временной зависимости. Описывается методика реконструкции математических моделей по рядам, которая иллюстрируются примерами из практики автора и его коллег.
Ключевые слова: Временные ряды, обработка сигналов, реконструкция динамических систем, эмпирическое моделирование, диагностика связей, приложения в биомедицине и науках о Земле.
Введение
Вид временного ряда - дискретной последовательности отсчетов - имеют сигналы на выходе практически всех современных измерительных устройств, использующих цифровую обработку, в том числе электроэнцефалографов, кардиографов, миографов и многих других медицинских приборов (пример {ц} на рис. 1). Форму ряда обычно имеют и записи статистических данных. С помощью сканирования в виде рядов может быть представлена информация о свойствах объектов, распределенных в пространстве. Ряды могут быть векторными и скалярными в зависимости от вида их элементов. Все имеющиеся подходы к обработке таких данных можно разделить на непосредственную обработку ряда и опосредованный путь с использованием дополнительной математической модели. Первый подход более традицио-нен - например, кардиолог непосредственно рассматривает кардиограмму, визуально оценивая форму зубцов; аналогично, непосредственно по записи сигнала (по временному ряду), только с помощью компьютера, можно оценить средние значения, спектры, корреляции или построить фазовые портреты, рассчитать меры их сложности и т.п. При втором подходе для анализа процесса используется его модель, описывающую динамику — изменение переменных во времени (в общем случае это - не модель структуры объекта или механизма его функционирования). Опосредованный
70
Рис. 1. Вид записей сигналов активного мозга, звукового, потенциала на коже, связанного с активностью мышц, давления в нефроне почки. Все они имеют вид временного ряда (см. справа увеличенный фрагмент графика и соответствующий численный ряд)
подход добавляет к достоинствам непосредственного анализа сигналов возможности прогноза поведения во времени или при изменении параметров, классификации сигналов и диагностики, хранения данных, оценки адекватности представлений об объекте или о механизмах его функционирования, измерения величин не поддающихся прямой регистрации и т.п.
О том, как создать по временному ряду такую математическую конструкцию, почему для сложных движений предпочтительны уравнения, а не функции времени, и какие специфические приемы можно использовать, изложено ниже.
Начнем с идеи и предыстории задачи реконструкции (раздел 1), затем рассмотрим типичные подходы и некоторые оригинальные приемы восстановления уравнений и некоторые поучительные примеры (раздел 2), а закончим иллюстрацией нескольких приложений описываемой техники на практике (раздел 3). Достаточно полный перечень публикаций по рассматриваемым вопросам потребовал бы слишком много места, поэтому предлагаем заинтересованному читателю такой список в доступной на сайте (www.nonlinmod.sgu.ru) книге Безручко Б.П., Смирнова Д.А. «Математическое моделирование и хаотические временные ряды» (Саратов: ГосУНЦ «Колледж», 2005), которая легла в основу обзора. Для иллюстраций используем результаты, ранее изложенные в статьях нашей научной группы из СГУ и СФ ИРЭ РАН; ряд уникальных экспериментальных данных для расчетов был представлен нашими партнерами: из Москвы И.Г. Моховым из Института физики атмосферы РАН (климатические данные) и из Германии П. Тассом из Института медицины Исследовательского центра города Юлих (медицинские данные).
71
1. Идея реконструкции по временному ряду
Предшественниками современных задач математического моделирования по экспериментальным рядам являются известные с середины XVIII века задачи аппроксимации - описания явной временной зависимости наблюдаемой величины П = f (t,c) некоторой функцией f, где с - P-мерный вектор параметров модели. Их можно интерпретировать как проведение кривой через экспериментальные точки на плоскости (t, п) или вблизи этих точек (рис. 2). Если качественный вид связи между величинами t и п - функция f (t, c) - известен с точностью до значений компанентов с, задача состоит в том, чтобы как можно точнее оценить параметры зависимости. В другой, более сложной, постановке задачи требуется найти функцию f, способную обеспечить прогноз поведения с минимальной погрешностью. В примере, представленном на рис. 2, напрашивается предположение о поиске f в виде перевернутой и смещенной параболы. Здесь действительно фигурирует парабола, поэтому три параметра модели в отсутствие шумов (рис. 2, а) точно определяются из системы алгебраических уравнений с известными левыми частями
П1 = Со + Citi + C2tf,
П2 = Со + Cit2 + С2Ц, (1)
П3 = Со + Cit3 + C2t^,
составленной для любых трех различающихся моментов времени ti, t2, t3. В стохастическом случае (рис. 2, б) модель ищется в виде
П = f (t,Co)+ Ш, (2)
где x - случайный процесс («шум») с нулевым средним1. Здесь при расчете параметров речь идет не о точном их значении, а о получении статистических оценок с, как можно более близких к со.
Рис. 2. Иллюстрация задачи - провести кривую заданного вида через экспериментальные точки (кружочки) на плоскости: а - шума нет, вид функции подобрать несложно; б - точки не ложатся в точности на простую кривую из-за случайной помехи, либо из-за сложного характера зависимости. Задача состоит в аппроксимации наблюдаемой зависимости простой функцией
1В статистике процедура создания модели по рядам известна под названием «идентификация систем».
72
Рис. 3. Сложная хаотическая зависимость (а) точно аппроксимируется квадратичным отображением (б)
С усложнением множества точек на плоскости решение становится более громоздким - для аппроксимации приходится использовать все более длинные ряды Пз = Со + cit + c2t2 + c3t3 + ... + cktk, но и эти возможности ограничены. Например, аппроксимация простой синусоиды десятком слагаемых обеспечивает хорошее приближение на интервале длиной всего в несколько периодов. Очевидно, что аппроксимация хаотического временного ряда, приведенного на рис. 3, а, требует несравненно еще более громоздкой конструкции. Но ситуация резко упрощается, если отказаться от аппроксимации ряда функцией времени, а искать уравнение, решением которого является рассматриваемый ряд2. Так, хаотический временной ряд (см. рис. 3, а) является решением очень простого нелинейного разностного уравнения
Хп+1 = ГХп(1 - Xn), (3)
где r - параметр, а n = 0,1, 2, 3,... - дискретное время. Это - известное точечное квадратичное отображение, на котором в 1970-е годы исследовались закономерности хаотической динамики, а его график - такая же, как в (1), - парабола (рис. 3, б), но на других по сравнению с рис. 2 осях (на рис. 3, б по осям отложены значения переменной в данный и последующий моменты времени). А раз так, то задача реконструкции модели в виде x(tn+i) = f (x(tn), c) на этом этапе не усложняется по сравнению с рассмотренным ранее примером аппроксимации явной временной зависимости (хотя добавляются трудности догадки - какие величины следует взять в качестве переменных).
При аппроксимации дифференциальным уравнением dx/dt = f (x, c) по осям откладываются динамические переменные и скорости их изменения. Так, для синусоиды, являющейся решением уравнения гармонического осциллятора, переменными являются x1 = X и x2 = dx/dt, а модельная система имеет вид
dx1
dt
dx2
dt
fi(xi,x2) = x2,
h(xi,x2) = —w2xi.
(4)
Имея временные ряды переменных xi ,x2 и скоростей их изменения dxi/dt, dx2/dt, можно провести подбор коэффициентов искомых функций.
2Уравнением называют запись задачи в виде равенства двух функций. Решение задачи означает нахождение зависимости переменной уравнения от времени, подстановка которой превращает равенство в тождество. Откровением XX века стало понимание способности системы трех и более одномерных нелинейных дифференциальных уравнений демонстрировать хаотическое поведение. Для дискретных систем это возможно уже при одной переменной.
73
Рис. 4. Схема процесса реконструкции уравнений по временным рядам
Итак, реконструкция уравнения проводится по следующей схеме (рис. 4): исходя из вида временных рядов и априорной информации выбирается тип уравнений (разностное или дифференциальное), их структура (вид функций, число переменных Xi, способ получения Xi из временных рядов наблюдаемых величин {п}), а затем по экспериментальным значениям, взятым из рядов переменных и скоростей, проводится подгонка коэффициентов любым из известных для этого методов. Чем больше нам известно об искомой модели, тем больше вероятность успеха. Наиболее проста ситуация, когда известно все, кроме значений коэффициентов в функциях, а наиболее сложна, если сведения об объекте моделирования отсутствуют (имеем «черный ящик»). Но на каждом из этапов, даже самом простом (правом на схеме) возникают свои специфические технические трудности, для преодоления которых разрабатываются специальные технические приемы - технологии (см. литературу на сайте www.nonlinmod.sgu.ru). Рассмотрим некоторые из них.
2. Некоторые технические приемы процедуры реконструкции
2.1. Специфические подходы к выбору структуры модели. Неопределенность в выборе структуры модели может быть уменьшена за счет использования априорной информации о свойствах системы или предварительного анализа рядов наблюдаемых. Проиллюстрируем это на двух важных примерах.
2.1.1. Восстановление характеристик систем с запаздыванием. Речь пойдет о системах, чья динамика описывается дифференциальными уравнениями с запаздывающими аргументами вида
eX(t) = —x(t) + F(x(t — то)), (5)
где скалярная переменная х(наблюдаемая) может содержать и шумы. В принципе, можно пытаться решать задачу, используя данные ряда х c помощью метода наименьших квадратов, минимизируя сумму Yln (X(tn) — ax(tn) — f (x(tn—т), c))2 ^min и рассматривая при этом время задержки т как дополнительный неизвестный параметр. Но существуют априорные предпосылки и технология реконструкции, опирающаяся на анализ расположения экстремумов во временных реализациях наблюдаемых колебаний. Статистический анализ временных интервалов, разделяющих экстремумы во временных реализациях различных модельных и реальных систем с запаздыванием, позволяет установить, что зависимость от величины т числа N пар экстремумов временной реализации, удаленных друг от друга на время т, имеет четкий минимум при времени, соответствующем времени запаздывания системы (рис. 5). Далее с помощью сделанной по этому минимуму оценке т ~ то можно найти оценку параметра инерционности е и аппроксимировать нелинейную функцию F.
74
Рис. 5. а - временная реализация уравнения Ике-ды (6) для £о = 1.0, ц0 = 20.0, то = 2.0: б - число пар экстремумов M (т), нормированное на общее число экстремумов ряда, Ммин(т) = M(2.0); в - восстановленная нелинейная функция. Численные эксперименты с добавлением измерительного шума показывают возможность реконструкции при уровнях шума до 20 % (по отношению среднеквадратичных амплитуд шума и сигнала)
Описанный статистический подход к определению времени задержки может быть распространен на системы с запаздыванием более высокого порядка, а также на системы с несколькими различными временами задержки. Он экономен по отношению к машинному времени и использует для восстановления нелинейной функции все точки временного ряда, что позволяет по коротким временным рядам восстанавливать нелинейную функцию F даже в случаях слаборазвитого хаоса. Статистический подход не критичен к присутствию шума в наблюдаемой. На рис. 5, в в качестве примера приведены результаты реконструкции по хаотической временной реализации (рис. 5, а) уравнения Икеды, описывающего динамику пассивного оптического резонатора,
ex(t) = —x(t) + ц sin(x(t — т0) — x0).
(6)
2.1.2. Реконструкция неавтономных систем: регулярное воздействие. Часто информацию о наличии регулярного воздействия можно получить из анализа спектра мощности наблюдаемого ряда, который для таких систем обнаруживает дискретные составляющие. В этом случае в уравнения модели включаются описывающие воздействие функции, явно зависящие от времени. Например, если воздействие гармоническое и аддитивное, целесообразна модель в виде
dD x dx dD 1x
~цр = f (*’p’-’ dtD~1 ’ C) + a C°s “ + b sin “t, (7)
где f - алгебраический многочлен.
В случае произвольного регулярного воздействия (сложного периодического и квазипериодического) удобно представление
dDx
div
f (x,
dx dt ’ ’
dD-1x
dtp-i, c) + g(tc),
(8)
где функция g(t) описывает воздействие и также содержит неизвестные параметры. Она может иметь вид тригонометрического многочлена
75
(9)
к Ki
j(‘) = EE cAj cos(jmit + ).
i=i j=i
Причем адекватные модели с тригонометрическими многочленами могут быть получены и при очень большом числе (сотни) используемых гармоник, тогда как при аппроксимации алгебраическими многочленами увеличение размера модели чревато неустойчивостью их траекторий. Эффективность подхода показана на примерах реконструкции уравнений по зашумленным хаотическим реализациям эталонных осцилляторов при различных видах воздействия: периодическом импульсном, периодическом с субгармониками, квазипериодическом. Кроме того, зависимость от времени можно вводить и во все коэффициенты многочлена f, чтобы описать неаддитивное воздействие.
2.2. Выбор переменных. Пусть проблемы получения экспериментального ряда решены полностью, и мы определились с типом модельного уравнения. Какой из вариантов выбора переменных предпочесть? Теоретически можно пробовать поочередно все мыслимые варианты, для каждого выбора аппроксимировать зависимость вида dx/dt = f (x, c) или xn+1 = f (xn, c) и оценить результат, но это не реально. Желательно было бы заранее выбрать относительно небольшое число наиболее подходящих вариантов и только для них проводить аппроксимацию. Для осуществления такого выбора предложен ряд процедур, опирающихся на очевидное требование: переменные динамической модели должны обеспечивать однозначность и непрерывность зависимостей значений стоящих в «левых частях» уравнений величин от значений динамических переменных. Обозначим «левую часть» модельных уравнений z и проведем дальнейшую часть процедуры аппроксимации рядов наблюдаемой {x(ti)} и «левых частей» (z(ti)}3. Далее нужно проверить, соответствуют ли близким значениям x близкие соответствующие значения z. Для этого область V, внутри которой содержится множество векторов {x(ti)}, разбивается на одинаковые «гиперкубические» ячейки со стороной 8 (рис. 6, а). Из них выбираются все ячейки si, , которые содержат более одного вектора x(ti), то есть имеют хотя бы
две точки в нижней клетке рисунка. Разность между максимальным и минимальным значениями z (одной из компонент вектора z) в пределах ячейки Sk назовем локальным разбросом ек. По величине максимального разброса емакс = maxi<k<M ек и графику емакс(8) оценивается пригодность величин x и z для глобального моделирования. Для построения глобальной модели переменные нужно выбирать так, чтобы график емакс(8) стремился к началу координат плавно, без изломов (жирная линия на рис. 6, б) для каждой из аппроксимируемых зависимостей Zk(x), к = 1,D.
При этом желательно, чтобы наклон графика емакс(8) был как можно меньшим, потому что для аппроксимации тогда достаточно использовать более простую модельную функцию, например, многочлен низкого порядка. На рис. 6, в, г это проиллюстрировано на простом примере - аппроксимации зависимости следующего значения наблюдаемой от предыдущего, когда наблюдаемая генерируется первой, второй или третьей итерацией квадратичного отображения x(tn+1) = к — x2(tn). График первой итерации «наименее осциллирующий», а потому наклон емакс(8) самый маленький. Именно в этом случае проще всего получить «хорошую» модель
3Для обыкновенных дифференциальных уравнений это делают численным дифференцированием ряда {х(^)}, а для отображений - сдвигом {x(fi)} на один шаг по времени.
76
штриховая линия - худший (неоднозначность или разрывность), тонкая линия соответствует зависимости с областями быстрого изменения; в - графики 1, 2 и 3-й итераций квадратичного отображения (помечены цифрами); г - графики £мжс(б) для зависимостей x(tn+i) от x(tn) в этих трех случаях
(действительно, достаточно многочлена второго порядка, тогда как для описания третьей итерации нужен многочлен восьмого порядка). Особенно различаются эти три случая по трудности реконструкции при наличие шумов.
Самый простой способ получения дополнительных переменных по ряду наблюдаемой это метод временных задержек. В случае скалярной наблюдаемой, согласно этому подходу, в качестве компонент вектора x(t) берутся ее последовательные значения, разделенные некоторым интервалом т (временем задержки): Xl(t) = n(t), X2(t) = n(t + т), ..., XD (t) = n(t + (D - 1)т).
Не менее популярен и метод последовательного дифференцирования, который состоит в использовании временных производных наблюдаемой в качестве динамических переменных: xi(t) = n(t), x2(t) = dn(t)/dt, ..., xD(t) = dD-1n(t)/dtD-1.
Его применение затруднительно при наличии измерительных шумов, так как шумовая «бахрома» имеет на графике n(t)большую крутизну, после дифференцирования уровень шума резко возрастает.
Имеющиеся шумы частично подавляются при использовании интегрирования для получения дополнительной переменной. Например, при использовании интегралов от наблюдаемой с переменным верхним пределом: x2(t) = J0 n(t)dt. Однако при интегрировании теряется часть информации о наблюдаемой.
77
С помощью взвешенного суммирования можно составить переменную из последовательных значений наблюдаемой
X2(t) = aon(t) + ain(t — At) + a2n(t — 2At) +
где ak - весовые коэффициенты. Возможны комбинации всех упомянутых подходов и другие способы получения рядов модельных переменных.
В заключение остановимся на очень популярном в настоящее время использовании в качестве переменой «фазы колебаний» наблюдаемой. Фазой гармонического сигнала x(t) = Acos(mt + ф0) называют аргумент косинуса ф = rnt + ф0. Существуют обобщения понятия фазы на более сложные сигналы, из которых наиболее распространен вариант, использующий построение аналитического сигнала. Аналитический сигнал z(t) = x(t) + iy(t) строится по исходному сигналу x(t) с помощью преобразования Гильберта
y(t) = p.v.
,<х x(t')dt' -<х n(t — t) ’
(10)
где P.V. означает главное значение несобственного интеграла. Это преобразование поворачивает фазы гармоник исходного сигнала на п/2. Фаза определяется как аргумент z(t), то есть как угол поворота радиуса-вектора на комплексной плоскости (x,y). Так, в случае гармонического сигнала x(t) = A cos(mt + ф0) получаем сопряженный сигнал y(t) = Asin(mt + фо). Фазовая траектория на плоскости (x,y) представляет собой окружность (рис. 7, а). Фазу определяют как угол поворота изображающего вектора относительно оси x. Для более сложного («похожего» на гармонический, но узкополосного) сигнала наблюдается вращение радиус-вектора z(t) не по окружности, а «почти» по окружности (рис. 7, б). Фаза нарастает в среднем со скоростью, равной средней угловой частоте колебаний. В случае сложных нерегулярных сигналов без ярко выраженного основного ритма может оказаться, что построение аналитического сигнала через преобразование Гильберта не обнаружит вращения около четко определенного центра (рис. 7, в). Тогда такое формальное определение фазы, бесполезно. Поэтому часто используют полосовую фильтрацию, чтобы сузить на сколько возможно полосу частот сигнала.
Рис. 7. а - траектория на комплексной плоскости для гармонического сигнала, его амплитуда и фаза; б -определение амплитуды и фазы негармонического узкополосного сигнала с помощью преобразования Гильберта; в - то же самое для широкополосного сигнала (фаза не является хорошо определенной величиной)
78
2.3. Оптимизация структуры модели. Одной из проблем реконструкции уравнений является чрезмерная громоздкость аппроксимирующих функций из-за использования длинных функциональных рядов, в которых многие члены лишние. Лишние слагаемые можно идентифицировать путем наблюдения за поведением оценок коэффициентов модели при реконструкции по различным участкам временного ряда, то есть по точкам, лежащим в различных областях фазового пространства. Параметры адекватной глобальной модели стационарной системы не должны зависеть от того, по каким участкам ее временного ряда ведется реконструкция. А коэффициенты при лишних слагаемых могут претерпевать, в зависимости от выбора участка, заметные изменения. Причем наиболее сильно это проявляется, если временной ряд содержит переходный процесс, так как при этом фазовой траекторией «осваиваются» существенно различные области фазового пространства.
Будем строить модель с многочленом высокого порядка K по последовательным сегментам временного ряда длиной W: {n(k-i)w+i, ■■■, n(k-1)w+w},
(k)
к = 1, 2, ...,L. Получим набор оценок коэффициентов c\j (рис. 8, а). Степень стабильности оценки каждого коэффициента j определим как модуль отношения эмпирического среднего {ci j) = (1/L)^^=1 с( к к стандартному отклонению
та сегмента реконструкции; б - изменение погрешности при исключении наиболее осциллирующих коэффициентов; погрешность минимальна при 20 исключенных слагаемых
79
y(1/L)YlU (Ck,j — {Ck))2. Исключим слагаемое, соответствующее наименее стабильному коэффициенту. Для упрощенной структуры модели повторим всю процедуру. Продолжая ее дальше, будем последовательно удалять «нестабильные слагаемые». Исключение слагаемых прекращается, когда качество модели перестает улучшаться. На рис. 8, б после исключения 20 слагаемых из модельного многочлена (K = 7) погрешность уменьшается на порядок по сравнению со стартовой структурой модели.
Перечень технических приемов можно было бы продолжить, но ограниченный объем статьи не позволяет сделать сколь либо подробный обзор.
3. Примеры приложений методов опосредованного анализа сигналов с помощью эмпирических моделей
Пожалуй, самым известным приложением является использование моделей для прогноза будущего по наблюдаемому временному ряду. У всех на слуху прогнозы погоды и климата, землетрясений, курсов акций и валют и пр. Положительный результат реконструкции модели (хорошее ее качество) может быть истолкован, как одно из доказательств правильности физических соображений, заложенных в ее структуру - оценка адекватности представлений, заложенных в модель. Это самодостаточная фундаментальная ценность, которая может привести впоследствии к разнообразным практическим приложениям.
Среди прочих полезных приложений можно выделить задачи выявления квазистационарных участков в нестационарном сигнале, прогноз бифуркаций слабонеавтономных систем, многоканальная конфиденциальная передача информации. Отметим также задачи восстановления сигнала внешнего воздействия по хаотическому ряду координаты неавтономной системы, классификации сигналов, управления техническими объектами, расчета характеристик аттракторов по коротким временным рядам, тестирования на нелинейность и детерминизм, проведение адаптивной фильтрации. Применительно к тем или иным профессиональным запросам могут оказаться наиболее полезными определенные «инструменты». В качестве примера рассмотрим следующие «инструменты», которые можно предложить биомедицине.
Анализ записей сигналов. Запись потенциалов мозга (электроэнцефалограммы - ЭЭГ), с мышц (миограммы), с сердца (кардиограммы), сделанные с электродов наряду с необходимой информацией содержит различные искажения - артефакты, мешающие выявлению патологий (рис. 9, а). Такие участки можно выделить с помощью предсказательной модели в виде D-мерных отображений x(tra+i) = f (x(tn), c). Построенная по участку с нормальной активностью такая модель предсказывает последующее поведение, а сделанная ошибка может быть использована как критерий отклонения от выбранной нормы (рис. 9, б). Наоборот, если эталонная модель составлена по участку с конфигурацией, присущей некоторой патологии, то малая ошибка прогноза будет свидетельствовать о наличии признаков такой патологии, а большая - об отличии от свойств заданного эталона.
Аналогичный подход, когда по различным участкам ряда реконструируются одинаковые уравнения, позволяет выделять интервалы стационарности - участки
80
Рис. 9. Электроэнцефалограмма с присущими ей артефактами (а) и графики погрешности, где выбросы соответствуют интервалам отклонения от нормы (б)
записи, на которых коэффициенты модели остаются неизменными. Сопоставление коэффициентов моделей, подогнанных на различных участках записи, позволяет различить интервалы стационарности в нестационарном процессе.
Набор коэффициентов базовой модели может использоваться для запоминания данных при диспансеризации. Такой подход позволяет ввести количественную меру изменения состояния здоровья за промежуток после предыдущих обследований, если в качестве меры изменения состояния использовать расстояние между значениями параметров модели в пространстве параметров (например, среднеквадратичное значение изменений всех параметров модели).
Анализ связей между элементами системы. Если структура адекватной математической модели ансамбля элементов известна и имеется конечный набор вариантов структуры связей, параметры которых следует установить, можно действовать простым перебором - искать варианты, которые обеспечивают наиболее точное воспроизведение наблюдаемой динамики.
В отсутствие априорной информации о структуре модельных уравнений
К. Грейнджером предложен следующий подход выявления причинно-следственных связей. Основная идея состоит в том, что по имеющимся временным рядам {x1(t1), ...,xi(tN)} и {x2(t 1), ...,x2(tN)} строятся прогностические модели - «индивидуальные» и «совместные». Существенное улучшение прогноза динамики первой системы при усложнении модели (учете значений переменной еще и от второго элемента) воспринимается как признак наличия воздействия второй системы на первую (если такого улучшения нельзя достичь усложнением индивидуальной модели).
Наиболее чувствительным и приспособленным к выявлению предпатологий является анализ фазовой динамики. Этот подход применим, если в исходных временных рядах имеют место ярко выраженные ритмы колебаний (в спектрах мощности имеется существенный подъем в узкой полосе частот). В этом случае можно недвусмысленно определить фазы ф1 и ф2 наблюдаемых процессов. Широкому классу колебаний адекватно представление в виде достаточно простой системы стохастических разностных уравнений
Ф1,2 (t + Т) - ф 1,2 (t) = f 1,2 (Ф 1,2 (t), Ф2,1 (t - Al,2)), (11)
где f1,2 - тригонометрические многочлены невысокого порядка; t - фиксированный
81
временной интервал, обычно равный меньшему из характерных периодов колебаний; Ai,2 - пробные времена запаздывания воздействий. Коэффициенты тригонометрических многочленов оценивают по временному ряду и по полученным оценкам рассчитывают силы воздействия систем друг на друга. Степень влияния с2^\ системы 1 на систему 2 определяется крутизной зависимости f от ф2, аналогично для с1^:
Недостатком рабочей формулы (12) является требование довольно большой длины исходного ряда (порядка 1000 условных периодов). Но это требование может быть смягчено на порядок, если ввести поправки п,2
зависящие от уровня шума, частот колебаний и длины временного ряда.
Перечисленные варианты оценок связанности полезны в приложении ко многим заболеваниям нервной системы, включая эпилепсию и болезнь Паркинсона, при которых имеет место патологическая синхронизация нейронов. Однако функциональная роль этой синхронизации при генерации паркинсоновского тремора (непроизвольных регулярных колебаний конечностей с частотами от 3 до 6 Гц) остается предметом обсуждений. Несмотря на неполную ясность механизма патологии, в качестве стандартной терапии при отсутствии эффекта от лекарств используется постоянная глубокая электрическая стимуляция мозга. Процедура небезобидна - не исключены побочные эффекты, а уменьшить меру воздействия можно лишь, увеличивая целенаправленность воздействия. В связи с этим актуальна задача определения характера связей между различными областями мозга и колебаниями конечностей пациентов.
Задача решалась нами путем анализа экспериментальных данных с использованием записей локальных потенциалов с глубинных электродов и сигналов акселерометров, закрепленных на конечностях (рис. 10). Информация об активности мозга была представлена записями локальных потенциалов от четырех глубинных электродов, вживленных в таламус или базальные ганглии. Данные из нескольких десятков интервалов спонтанного тремора для каждого пациента были получены на факультете стереотаксической и функциональной нейрохирургии Университета Кельна и в Институте нейронауки и биофизики Исследовательского центра Юлиха.
С помощью нелинейного метода, основанного на моделировании фазовой динамики, была выявлена двунаправленная связь между процессами. Оказалось, что воздействие субталамического ядра на конечности является нелинейным эффектом и демонстрирует запаздывание в 200-400 мс. Обратное воздействие (конечности -субталамическое ядро) носит линейный характер с запаздыванием не более 50 мс. Такое взаимодействие обнаружено по экспериментальным данным впервые и означает, что субталамическое ядро является звеном «кольца обратной связи», определяющим колебания конечностей, а не просто пассивным приемником сигналов.
(12)
Y1 —>2,2—> 1 — с1—2,2—1 + r 1,2
,2
82
Рис. 10. Записи: сигнала с акселерометра (а), локального потенциала на электроде в мозге (б) и соответствующих им спектров мощности S (в). Оценки связанности на фоне доверительного интервала (д) при болезни Паркинсона: мозг ^ рука (г), рука ^ мозг (е)
Другие примеры. Методы обработки рядов, опирающиеся на моделирование динамики, позволяют рассматривать данные различной природы. В качестве примера в заключение рассмотрим результаты оценок связанности (по Грейнджеру) для анализа влияния солнечной активности на глобальную приповерхностную температуру (ГПТ) Земли, сделанных И.Г. Моховым и Д.А. Смирновым. Использовались две версии сигналов, отражающих солнечную активность - по Хойту и Лину4. Установлено наличие влияния Солнца на ГПТ без запаздывания или с запаздыванием в пределах 1 года. При этом в первой трети XX века солнечная активность определяет 5% дисперсии процесса ГПТ, а во второй половине XX века - 25 %. В короткий период 1925-1945 годов вполне могло происходить резкое изменение каких-либо характеристик процессов. Возможно, что с середины XX века воздействие Солнца усиливается до сих пор (вывод только по Хойту). Можно предположить, что характерный масштаб вариаций характеристик воздействия Солнца на ГПТ составляет около 30-50 лет. Все результаты анализа более стабильны при анализе по данным, полученным по методике Хойта.
4Hoyt D.V., Schatten K.H. The Role of the Sun in Climate Change. Oxford: Oxford Univ. Press, 1997. 279 pp.
Lean J., Rottman G., Harder J., Kopp G. Source contributions to new understanding of global change and solar variability // Solar Physics. 2005. Vol. 230. P. 27.
83
Заключение
Опосредованный подход к обработке сигнала опирается на математическую модель. Моделирование, по-видимому, всегда останется в значительной степени искусством, но могут быть выделены некоторые общие принципы и частные рецепты, позволяющие повысить шансы на получение «хорошей» модели. Ряд таких соображений, изложен в данной лекции. Но сказанное является лишь «экскурсом в...», а не исчерпывающим изложением техники эмпирического моделирования.
Выражаю искреннюю благодарность коллегам из СГУ и СФ ИРЭ РАН, чьи результаты использовались в лекции: Диканеву Т.В., Караваеву А.С., Селезневу Е.П., Сидак Е.В., Смирнову Д.А., Сысоеву И.В., Пономаренко В.И., Прохорову М.Д.
Изложенные результаты получены в рамках исследований, которые проводились при финансовой поддержке гранта РФФИ (08-02-00081) и Программы Президиума РАН «Фундаментальные науки - медицине».
Поступила в редакцию 4.08.2009
MODELING FROM TIME SERIES AND APPLICATIONS TO PROCESSING OF COMPLEX SIGNALS
B.P. Bezruchko
Signals obtained from most of real-world systems, especially from living organisms, are irregular, often chaotic, non-stationary, and noise-corrupted. Since modern measuring devices usually realize digital processing of information, recordings of the signals take the form of a discrete sequence of samples (a time series). The present paper gives a brief overview of the possibilities of such experimental data processing based on reconstruction and usage of a predictive empirical model of a time realization under study. The technique of reconstruction of mathematical models from time series is described and possibilities of the approach are illustrated with examples from the author's and his colleagues' experience.
Keywords: Time series, signal processing, reconstruction of dynamical systems, coupling diagnostics, empirical modelling, biomedical and geophysical application.
Безручко Борис Петрович - родился в 1946 году. Окончил физический факультет Саратовского госуниверситета (1969). Доктор физико-математических наук (1995). Заведующий кафедрой биомедицинской инженерии и динамического моделирования Саратовского госуниверситета, заведующий лабораторией моделирования в нелинейной динамике СФ ИРЭ РАН. Опубликовал более 100 статей в научных журналах и 2 монографии (в соавторстве). Область научных интересов: радиофизика и электроника, нелинейная динамика, моделирование по временным рядам с приложением к задачам физиологии и медицинской диагностики, физический эксперимент.
410012 Саратов, ул. Астраханская, 83
Саратовский государственный университет им. Н.Г. Чернышевского E-mail: [email protected]; [email protected]
84
