CC BY
83
Об одном методе моделирования нестационарных динамических систем и процессов
Ю.П. Драница
Судоводительский факультет МГТУ, кафедра высшей математики и программного обеспечения ЭВМ
Аннотация. Разработан корректный метод оценки нестационарности многомерных динамических процессов и ее учета при построении линейных математических моделей. Теория метода основана на аппроксимации процесса системой линейных разностных уравнений в прямом и обратном времени, при этом не требуется знания ковариационной матрицы. Метод является обобщением теории Дж. Бурга, разработанной для одномерного и стационарного случая. Выполнена оптимизация алгоритма за счет использования свойств матриц Теплица и рекурсии Левинсона.
Abstract. The correct method of an estimation of nonstationarity of the multidimensional dynamic processes and its taking into account when constructing linear mathematical models has been worked out. The theory of the method is based on approximation of the process by the linear differential equations system in direct and inverse time, and the knowledge of a covariance matrix is not required. The method is the generalization of the J. Burg theory, worked out for the one-dimensional and stationary case. The optimization of the algorithm has been executed by using properties of Teplitz matrixes and Levinson recursion.
1. Введение
Математическое моделирование геофизических систем и процессов является важным инструментом научных исследований. Однако на практике результаты моделирования иногда оказываются совершенно неприемлемыми, поэтому формальное использование этого метода не гарантирует получения достоверных решений. Требуется, чтобы принципы моделирования достаточно точно соответствовали модели данных, которая, в свою очередь, должна наиболее полно отражать физическую сущность явления.
Имеется несколько подходов к решению задачи улучшения качества математической обработки. Один из них состоит в вовлечении в обработку не только результатов непосредственных измерений, но и всей имеющейся априорной информации, которая позволяет в той или иной мере нейтрализовать недостатки принятой математической модели. Другой подход, не исключающий, а дополняющий первый, состоит в разработке все более точных и полных математических моделей, которые позволяют сблизить формальную математику с физикой решаемой задачи.
Известно, что большинство геофизических систем и процессов обладают ярко выраженной нестационарностью, или трендом. Однако в настоящее время не существует теорий или методов, в достаточной мере оценивающих и учитывающих это явление. Применение в этом случае хорошо разработанного стационарного приближения приводит к большим ошибкам в решении, а иногда даже к его расхождению. Поэтому задача корректного моделирования нестационарных систем и процессов является актуальной. Одному из возможных методов решения этой проблемы и посвящена настоящая работа.
2. Современное состояние проблемы 2.1. Исходные положения
Пусть наблюдения за некоторой динамической системой или процессом проводятся в дискретные моменты времени ti, i = 1,2,...,N, где N - общее число наблюдений. Измерения за состоянием системы характеризуется вектором S(ti) c размерностью m. Следовательно, совокупность наблюдений образует матрицу S размерностью mxN. Отметим, что вместо независимой переменной типа времени могут использоваться другие упорядоченные параметры, например, глубины, высоты и т.п., однако без ограничения общности будем считать, что независимой переменной является время.
Обычно при измерениях, дискретизации и квантовании данных возникают различного рода ошибки. В этом предположении процесс измерений может быть представлен следующей моделью
S(ti) = U(t)+E(t), (1)
где Щ((,) - вектор истинных значений параметров состояния системы в момент времени Е(() -суммарный вектор ошибок измерений в тот же момент времени; размерность вектора Е(^) равна т.
Выражение (1) представляет собой линейную модель процесса наблюдения и указывает на независимость сигнальной части измерений Щ(/,) и помехи Е(^). Сигнальная часть измерений, в зависимости от физической постановки задачи, может считаться либо неизвестным детерминированным вектором, либо набором случайных величин. Помехи, с одной стороны, представляют ошибки измерений, а с другой - отражают неполноту и ограниченность модели данных (1). В общем случае под помехами будем понимать любую информацию, которая не является полезной для отыскания вектора истинных значений параметров состояния системы в рамках принятой модели. Иногда помехи называют "шумом" наблюдений. Неполнота и ограниченность модели (1) может, в частности, заключаться в том, что реальные данные S(t1) содержат также систематические компоненты, или тренды параметров. В этом случае более адекватной моделью наблюдений будет следующая
S(t) = Л«д + У(Г)+Е(Г), (2)
где Л^), У(^) - соответственно вектора трендовой и "сигнальной" составляющих данных размерности т; Е(/,) - по-прежнему вектор шума.
Относительно суммарного вектора ошибок Е(^) принято считать, что он образует квазислучайную (коррелированную) совокупность данных. Трендовая составляющая отражает временную нестационарность процесса измерений и может иметь достаточно сложную природу. Например, она может состоять из систематического ухода показаний измерительных приборов, связанных с техническими особенностями этих устройств (приборные ошибки), и смещений, которые отражают реальный тренд наблюдаемого объекта. В этом случае вместо (2) должна быть использована более сложная модель, учитывающая данную реальность. Будем, однако, считать, что все приборные ошибки не имеют тренда и сосредоточены в "шумовой" компоненте сигнала. Это не умаляет общности рассуждений, так как при необходимости всегда может быть применена более адекватная модель наблюдений.
2.2. Основные проблемы моделирования
Из приведенных выше рассуждений видно, что разделение наблюденных данных на компоненты носит достаточно условный характер. Однако с методической точки зрения такое разделение является плодотворным, отражающим реальную дифференциацию данных. Особенно это эффективно, когда разделение может быть проделано на основе каких-либо априорных соображений, или исходя из физики процесса. Однако в реальных условиях такая информация, как правило, отсутствует, и реализацию моделей (1-2) приходится выполнять непосредственно по наблюденным данным. Рассмотрим основные особенности постановки задачи определения и учета параметров тренда при моделировании динамических систем.
Допустим, что по наблюдениям типа (1) или (2) требуется построить математическую модель динамической системы. Наиболее простой путь решения этой задачи заключается в использовании линейного метода. Математической основой такого подхода являются хорошо развитая теория линейных стационарных операторов и преобразование Фурье, которые позволяют ввести понятия передаточной функции, спектра, эргодичности и пр. Известно, что исчерпывающей информацией о линейной динамической системе является матрица ковариаций вектора состояний, позволяющая синтезировать ее математическую модель, которая ясно отражает такие ее физические свойства, как распределение энергии по частотам, устойчивость, наличие периодичностей, характер их затухания или возбуждения и другие параметры. Следовательно, основной задачей при линейном моделировании является корректная оценка ковариационной матрицы вектора состояний, в том числе его трендовой, "сигнальной" и "шумовой" компонент.
Линейные модели обычно подразделяются на модели с постоянными коэффициентами и с коэффициентами, зависящими от времени. Первая группа предназначена для описания стационарных процессов, а вторая - нестационарных. Анализ литературных данных показывает, что полная и законченная теория в настоящее время существует для стационарного приближения. Попытка учета нестационарности сопровождается резким усложнением теории метода и, соответственно, вычислительных алгоритмов. Причина этого явления заключается, вероятно, в том, что учет нестационарности, по сути дела, эквивалентен переходу к нелинейному моделированию, так как при этом в той или иной степени затруднено использование основных постулатов линейного подхода
(например, нельзя применять принцип эргодичности). Исключение составляют случаи, когда параметры нестационарности могут быть определены из физики процесса и аппроксимированы достаточно простыми аналитическими выражениями. Однако в практических приложениях такие ситуации встречаются исключительно редко.
Следует также отметить, что понятия "стационарность" или "нестационарность" - это достаточно абстрактные философские категории, сущность которых может проявляться только в статистическом смысле и для бесконечных выборок данных. Следовательно, для конечного набора данных не существует каких-либо объективных критериев, позволяющих установить стационарность либо нестационарность анализируемого процесса, и имеет смысл говорить только лишь о их некотором противопоставлении. Необходимо более конкретное наполнение этих понятий. На практике это означает, что по ограниченным по длине выборкам данных локальную изменчивость следует определять, исходя из доступного идентификации диапазона частот процесса, а более низкие частоты следует считать трендом. Тогда процесс можно считать стационарным, если энергия этих низкочастотных составляющих намного меньше энергии идентифицируемых частот.
Рассмотрим, к чему может приводить неучет тренда.
Эта задача тесно связана с проблемой влияния скрытых систематических ошибок на параметры идентифицируемой модели. Пусть полная модель некоторой совокупности данных описывается, например, выражением (2). Будем считать, что часть этой модели нам неизвестна и представляет собой скрытую систематическую ошибку. Это означает, что вместо соотношения (2) мы используем неполную модель (1), где ¿7(4) = А (4) + К(4). В работе (Макузе, 1989) показано, что скрытые ошибки приводят к систематическому смещению оценок параметров модели относительно истинных значений. Практические следствия этого могут быть весьма существенными, способными создать совершенно неадекватные представления о свойствах и параметрах изучаемого явления. Таким образом, задача объективного выявления нестационарности и ее учета является актуальной как с точки зрения разработки моделей с переменными коэффициентами, так и с постоянными для выработки критерия адекватности модели натурным данным.
Несмотря на всю важность влияния нестационарности на параметры модели динамического процесса, вопросам ее оценки и учета в литературе уделяется незаслуженно мало внимания. Зачастую данная проблема или вообще игнорируется, или рассматривается вскользь. При теоретических исследованиях задач линейного моделирования чаще всего исходят из предположения о стационарности процессов, или возможности приведения их к этому виду какими-либо другими, независимыми от разрабатываемой теории способами. При практической реализации алгоритмов линейного моделирования, при явном тренде данных, обычно рассматривается тренд процесса, в основном, по математическому ожиданию. Учет же тренда по дисперсии и ковариации, в силу объективных причин, обычно не выполняется. Отсутствуют также какие-либо оценки погрешностей, обусловленных принятыми приближениями. Вполне очевидно, что для построения моделей с высокой точностью необходимо сближение формального математического аппарата с физикой изучаемого явления, и учет тренда необходим. Последнее может быть достигнуто разработкой более полных и точных математических моделей.
Между тем эта проблема актуальна еще и потому, что оценка тренда представляет большой самостоятельный интерес, по крайней мере, в двух аспектах. Во-первых, тренд позволяет в той или иной мере судить о низкочастотной части спектра процесса, выявление которой в силу малости выборки объективно невозможно. Во-вторых, эти оценки совершенно необходимы при экстраполяции на период времени, где влияние нестационарности существенно, и ее учет объективно необходим.
2.3. Современные методы оценки тренда
Рассмотрим способы оценки тренда параметров системы классическими методами. Из приведенных выше рассуждений следует, что оценка нестационарности может быть сведена к оценке тренда матрицы ковариаций вектора состояния. В общем случае необходимо рассматривать тренд по математическому ожиданию, по дисперсии и по ковариации. И в связи с этим практическое решение этой задачи может быть разбито на ряд этапов. На первом шаге оценивается тренд процесса по математическому ожиданию. Для этого обычно прибегают к аппроксимации временного ряда степенными или тригонометрическими полиномами невысокого порядка. На основе этого приближения исходная информация на участке анализа центрируется, обеспечивая таким образом стационарность данных по математическому ожиданию. Если используется прогноз за пределы участка анализа, то тренд по математическому ожиданию экстраполируется также на участок
прогноза. После этого можно оценить дисперсию и ковариацию на различных отрезках временной оси, а затем аппроксимировать наблюдаемый тренд теми или иными аналитическими выражениями.
Основная возникающая при этом проблема заключается в том, что для того, чтобы сделать эти оценки с достаточной для практических нужд точностью, требуются большие объемы информации. Так, например, в (Котюк и др., 1967) показано, что в случае разложения функции автокорреляции (ФАК) по системе затухающих косинусоид ошибки ее оценивания представляются суммой ошибок по каждой гармонике. При этом надежное оценивание ординат ФАК возможно лишь для тех составляющих процесса, период которых в 5-10 раз меньше продолжительности наблюдений. Кроме того, ошибка оценивания возрастает с увеличением временного запаздывания между функциями, следовательно, достоверные оценки ФАК могут быть получены только при небольших сдвигах. Это означает, что по ограниченным по длине данным нельзя с достаточной для практических целей точностью оценить не только тренд процесса, но и его наиболее низкочастотные компоненты. В то же время синтез модели только по достоверным оценкам ковариационной функции неизбежно приводит к ее недостаточному разрешению в области низких частот, что снижает ценность получаемых решений. Однако построение линейной динамической модели возможно и без знания ковариационной функции. Рассмотрим эту задачу более подробно.
3. Предлагаемый метод решения задачи
Примем к рассмотрению некоторую дискретную динамическую систему (или процесс), состояние которой меняется во времени. Для линейного случая поведение такой системы можно описать следующим стохастическим разностным уравнением
Х(4) = ЪЛы(Гк) Х(Гк-) + /(4), ■ = 1,2,...,/ , (3)
где Х(4) - т-мерный вектор состояния системы (фазовые переменные); Лы(4) - заданные матрицы размерности тхт (переходные матрицы);- неизвестный т-мерный возмущающий вектор.
В общем случае фазовый вектор X определяется как множество таких переменных, которые вместе с известными входными сигналами необходимы и достаточны для нахождения отклика системы. Матрицы Л являются матрицами частных производных вида 5Х(4)/5Х(4_¿) ■ = 1,2, ... ,/, которые, возможно, зависят от времени. Возмущающий вектор /(4) представляет, с одной стороны, воздействие на систему внешней среды, а с другой - ошибки, связанные с неадекватностью модели натуре.
Уравнение (3) называется дискретной моделью состояния системы. Это линейное разностное уравнение порядка тх/. С другой стороны, эта модель является уравнением многомерной линейной авторегрессии, или одношаговым линейным предсказателем. Это уравнение позволяет предсказать значение фазового вектора Х(4) в момент времени 4 по его значениям в прошлые моменты времени к I, ... , tk_/. Число участвующих в прогнозе прошлых значений фазового вектора называется лагом модели. Модель (3) представляет так называемый физически реализуемый линейный многомерный цифровой фильтр. Его мерность определяется числом компонент фазового вектора.
В практических приложениях иногда рассматривают физически нереализуемые фильтры, которые позволяют предсказать текущее значения вектора по его / будущим значениям. Модель такого фильтра может быть представлена выражением, аналогичным соотношению (3)
Х(4) = ЕгА,(4) Х(4+) + /(4), ■ = 1,2,.,/. (4)
В частном случае, когда матрицы Лы, Лр■ не меняются во времени, модели (3-4) описывают стационарную во времени систему. При одномерном фазовом векторе (т = 1), модели (3-4) упрощаются. Векторы Х(4) превращаются в скалярные величины Х(^), вектора/(4) - в случайные функции времени Д4), а совокупность матриц Л - в вектор строки. Модель (3) часто называют перспективным предсказанием, а (4) - ретроспективным.
Задача идентификации моделей (3-4) заключается в определении матриц Л¿(4) и лага / по априорным или экспериментальным данным. При известной матрице ковариаций вектора Х, решение этой задачи хорошо известно (Клаербоут, 1981). Если фазовый вектор непосредственно измеряем, то его матрица автоковариаций может быть оценена и по экспериментальным данным. Но, как показано выше, существенным ограничением такого подхода является то, что достигнутая при этом точность является достаточно высокой только при небольших лагах модели. Однако практически лаг модели заранее, как правило, неизвестен, и должен быть оценен по экспериментальным данным. Все эти
проблемы сильно ограничивают возможности определения параметров моделей (3-4) по методу непосредственной эмпирической оценки ковариационной матрицы.
Для одномерного случая и в стационарном приближении синтез моделей (3-4) возможен и при неизвестной автоковариационной матрице (Сильвия, Робинсон, 1983; Клаербоут, 1981). Основная идея этого метода была разработана Дж. Бургом для расчета спектров высокого разрешения и заключалась в симметричном использовании как перспективных, так и ретроспективных предсказаний. Важным достоинством этого подхода, по сравнению с классическими схемами, является возможность построения моделей с лагами, превышающими область достоверных оценок ФАК, рассчитанных по общепринятым методам. Последнее достигается благодаря статистически более достоверной экстраполяции параметров модели. Преимущества этого метода особенно полезны при работе с короткими выборками, на которых общепринятые статистические методы часто оказываются неприемлемыми геофизически. Эта ситуация хорошо была изучена Дж. Бургом, в результате чего он использовал некоторые специфические свойства матриц Теплица для разработки своей процедуры максимально-энтропийного оценивания спектра и параметров одномерной линейной и стационарной модели. Вычислительная эффективность метода основана на рекурсии Левинсона, которая позволяет строить очень экономичные алгоритмы. Отметим, что рекурсия Левинсона для временной области имеет такое же важное значение, как и быстрое преобразование Фурье для частотной. Полученная в ходе вычислений информация может быть использована затем для более корректной эмпирической оценки корреляционной матрицы фазового вектора.
Распространим идею Бурга на синтез моделей (3)-(4) по экспериментальным данным на многомерный случай. При этом предполагается, что ковариационная функция неизвестна, а моделируемый процесс нестационарен. Синтез модели будем осуществлять методом последовательных итераций, с увеличением на каждом шаге лага модели на единицу. Рассмотрим перспективное предсказание (экстраполяцию вперед) и ретроспекцию (экстраполяцию назад). Эти экстраполяции будем помечать символом "€". Согласно рекурсивному методу Теплица, сначала будем находить уравнения авторегрессии (3)-(4) единичного лага, затем, используя найденные соотношения, строим модель с задержкой на 2 единицы, и т.д. пока не будет достигнута желаемая точность аппроксимации, или получаемые приближения уже не будут дальше улучшаться. Достигнутый при этом шаг итерации примем за эффективный лаг модели. Для реализации этой идеи на каждом шаге итерации будем вычислять достигнутую погрешность прогноза. В заключение отметим также, что рекурсия Теплица является классической рекурсией, встречающейся в теории полиномов, ортогональных на единичном круге. Итак, на /-ом итерационном шаге имеем
х€Ы(4) = Хр,(4-1+1) Сь/ + Bi(t|c), Х€рг(Ь-г+1) = Хы^к) Ср/ + Рг(Ь-г+д. (5)
Здесь Сы, С^ - матрицы частных производных размерности шхш (переходные матрицы); значки "Ь", "р" обозначают принадлежность матрицы соответственно для перспективного и ретроспективного прогноза; второй нижний индекс у этих матриц означает номер итерации; В,(4), Р,(4) - вектора квазислучайных невязок с матожиданием Е(В,) = Е(Рг) = 0 и размера т; к = /,(г'+1),...,^ / = 2,3,...,/ - номер итерации; 0 -нулевой вектор размера т.
Из формулы (5) видно, что число отсчетов между прогнозируемой величиной и предиктором равно |/-1|, поэтому начальная итерация имеет номер /„ = 2. Определим начальные условия (при / = 1) следующим образом: СЬ1 = Ср1 = Е, где Е - единичная матрица размера шхш. Выражение (5) является одношаговым предсказателем, использующим для прогноза текущего значения вектора одно его прошлое (будущее) значение, сдвинутое на ((-1) временных тактов. Представим соотношения (5) в матричном виде
Х€Ы = Хр гСЬ г + В /
Х€р / = ХЫСр / + Р,, (6)
где Хы, Хр,, Вг, Р, - матрицы размера (Ж-/)хт.
Проблема состоит в определении оптимальных (по некоторому критерию) значений матриц Сы, Срг. При (N-1) < т проблему (6) можно рассматривать как типичную задачу метода наименьших квадратов (МНК) с критерием оптимальности в виде суммы квадратов невязок.
J1 = min{(Xw - Xpl Cbi)' (Xb, - Xpi Cw)},
J2 = min{(Xp, - Хы Cp) (Xp, - Хы Cp,)}, (7)
где символом "' " обозначена операция транспонирования.
Нам необходимо найти такие оценки Cbi, Cpi параметров модели (6), которые удовлетворяли бы этому условию. Для раскрытия экстремума функций J1(Cbi), J2(Cp,) воспользуемся правилами векторного дифференцирования (Pao, 1968). Пусть Z, A, B - матрицы или вектора. Тогда имеют место следующие соотношения:
d(B'Z)/ dZ = d(ZB)/ dZ = B, d(ZZ)/ dZ = 2Z, d(B'AZ)/ dZ = A'B, d(Z'AZ)/ dZ = 2AZ. (8)
Используя эти соотношения, находим необходимое условие экстремума
5(J1)/ 8Z = - Xb Xp, + X\ Xp Cb, = 0, 5(J2)/ dZ = - Xv Xb , + X'b Xb Cp, = 0,
откуда имеем
Xp Xp Cbi = XbiXpi, Xb Xb iCp i = Xp Xb,. (9)
Условием минимума квадратичных форм J1, J2 означает, что матрицы
52(J1)/5Z2 = XpXp„ 32(J2)/ 5Z2 = Xb Xb,
должны быть положительно определены.
Выражение (9) представляет собой матричную запись так называемой нормальной системы уравнений. Если матрицы Xbi, Xpi неособенные, то системы (9) имеют единственные решения
Cbi = (XpiXpi) X'biXpi, Cpi = (XbiXbi) XpiXbh (10)
где символом "-1" обозначена операция вычисления обратной матрицы.
Множители правой части соотношения (10) представляют матрицы дисперсий и взаимных ковариаций между будущими и прошлыми значениями вектора состояний. Как следует из формулы (5), дисперсии, вычисляемые по предыдущему выражению, рассчитываются на разных временных базах (для перспективного прогноза используется база t¡ < tk < tN.i+1, а для ретроспективного - ti+1 < ti < tN). Временная несовмещенность этих баз зависит от номера итерации и в целом составляет 2i отсчетов данных. Взаимные ковариации, вычисляемые для различных временных задержек, также пропорциональны номеру текущей итерации. Следовательно, при увеличении лага модели уменьшается число используемых в вычислениях данных.
Упростим выражение (10). Для этого введем следующие обозначения, смысл которых легко установить при анализе предыдущих формул:
XbiXpi = Rbph XpiXbi = Rpbb XpiXpi Dpi,
XbXbi = Dbl. (11)
Тогда соотношения (10), с учетом обозначений (11), могут быть переписаны в более компактном виде
Cbi = (Dp¡)) 1 Rbpu
Cpl = (Dbi)-1 Rpbi. (12)
Отметим, что 5-ая строка матриц Cbi, Cpi является строкой весовых коэффициентов для прогноза 5-ой составляющей текущего фазового вектора по информации, задержанной на (i-1) временных тактов.
Известно, что решение (12) может иметь вычислительные трудности при обращении матриц Б, если они вырождены или плохо обусловлены. Это является достаточно серьезным ограничением, которое может разрушить сходимость итерационного процесса. В работе (Химмерблау, 1975) изложен способ, позволяющий улучшить обусловленность матрицы Б, в котором предлагается обращать матрицы вида (Б + ЛЕ), где Е - единичная матрица, Л — параметр регуляризации. Известна также модификация приведенного выше метода, в котором обращается следующая матрица (Б + ЛБ^), где Б^ -диагональная матрица, составленная из соответствующих элементов матрицы Б. Величина параметра Л находится либо путем оптимизации его значения, либо задается априорно. В последнем случае параметр регуляризации Л обычно принимается равным 0.01^0.05, т.е. порядка нескольких процентов. Другим способом решения этой проблемы является сингулярный анализ (Троян, Соколов, 1989).
Имея оценки переходных матриц (12), можно выполнить прогноз фазового вектора по следующим формулам:
где первый индекс у векторов В, Р означает дискретное время, а второй - номер итерации. В результате получим ошибки предсказания и ретроспекции, которые следуют из системы уравнений (5):
при начальных условиях Вк1 = Х(4), Рк1 = Х(4), к = 1,2,...,N.
Поясним физическую сущность введенных выше соотношений. С учетом формулы (11) и невязок (14), матрицы Бы, Бр/ интерпретируются как ковариационные матрицы ошибок предсказания и ретроспекции (с размерностями тхт). Вектора В, и Р, интерпретируются, соответственно, как ошибки предсказания и ретроспекции на /-ом шаге после учета линейных связей предыдущих (/-1) итераций. Поэтому матрицы ЯЬр/ и ЯрЫ представляют взаимные функции частных автоковариаций составляющих фазового вектора (ВФЧАК) предсказания и ретроспекции на (/-1)-ом сдвиге. И в этом случае, как видно из формул (12), матрицы Сы, Ср/ можно трактовать как нормированные значения ВФЧАК. Мы используем понятие частных автоковариаций, чтобы подчеркнуть тот факт, что на /-ой временной задержке учитываются только те линейные связи, которые еще остались после устранения линейности для более ранних запаздываний.
Анализ используемых соотношений показывает, что между перспективным и ретроспективным прогнозом имеются определенные функциональные связи, которые обеспечивают устойчивость и сходимость алгоритма. Такая взаимообусловленность моделей прямого и обратного времени обеспечивается своего рода обратными связями, которые имеются в формулах (13-14). Под обратными связями здесь понимается то, что при вычислении невязки перспективного (ретроспективного) прогноза используются ошибки как для прямого, так и для обратного времени, полученные на предыдущей итерации. Для реализации этой идеи алгоритм хранит соответствующие невязки в двух относительно независимых матрицах. С одной стороны, обратные связи обеспечивают монотонную сходимость итераций и согласование моделей между собой. С другой стороны, гарантируется минимальное отличие моделей перспективного и ретроспективного прогноза. Для стационарного случая теоретически эти две модели должны быть эквивалентными, а практически - должны иметь только статистически незначимые отличия. В этом случае обнаружение больших отклонений между ними сигнализирует о нестационарности параметров процесса.
Важно отметить, что при построении модели алгоритм не использует никакой априорной информации, а только текущие данные, которые содержит временной ряд, т.е. метод не использует статистически необоснованной информации. Более того, оптимизация производится на каждом итерационном шаге, поэтому итоговую модель можно рассматривать как кусочно-линейную. Такая методика позволяет учитывать конкретную статистическую структуру обрабатываемой информации на каждой временной задержке, что способствует построению более гибких и точных моделей. Раздельное построение моделей для прямого и обратного времен также способствует более дифференцированному учету текущей информации. В то же время классический МНК использует уже готовые оценки дисперсионных и ковариационных матриц, которые, как показано выше, могут
Вк/ = Р(к-/+1)(/-1) СЬЬ Р(Ы+Г)1 = Bk(i-11)Cpi,
(13)
Вк/ = ВЙ>1) _ Р(к-/+1) (¿-1) CЬ¿, Р(к-/+1) / = Р(к-1+1)( /-1) - Вк(/-1) Ср/
(14)
содержать значимые ошибки. Кроме того, оптимизация по классическому МНК выполняется только один раз, т.е. обеспечивает средневзвешенный оптимум по всему набору параметров модели и временным сдвигам. Разная статистическая структура прямого и обратного времени здесь также не учитывается. Таким образом, разработанная методика позволяет более тонко учитывать структуру данных и отражать ее в формируемой модели. В этом, по нашему мнению, заключаются принципиальные отличия разработанной методики от классического подхода.
Для стационарного процесса теоретические дисперсии являются постоянными, а между ВФЧАК существуют простые зависимости (Клаербоут, 1981). В целом все эти условия можно описать следующими соотношениями:
Db = Dp, (15')
Rbp = R'pb, (15")
Cb = C V (15'")
Рассмотрим одномерный вариант метода. В этом случае матрицы из выражений (15'-15"') вырождаются в числа, при этом ВФЧАК трансформируется в ФАК. Известно, что ФАК в одномерном случае является четной функцией, и, следовательно, для любого процесса (с трендом или без него) выполняется равенство Rb = Rp. При практической реализации алгоритмов в классическом варианте выполнение соотношения (15') достигается за счет усреднения дисперсий ошибок предсказания и ретроспекции. Анализ формулы (12) показывает, что при выполнении равенств (15'-15"), автоматически выполняется и соотношение (15"'). Многомерный стационарный случай не имеет в этом отношении принципиального отличия от одномерного, следовательно, и для него выполняются условия (15'-15"').
Рассмотрим теперь нестационарный процесс сначала в одномерном варианте. В этом случае теоретические дисперсии уже являются функциями времени и, следовательно, равенство (15') должно нарушаться, т.е. теоретические дисперсии предсказания и ретроспекции уже не могут считаться равными. Практические оценки этих дисперсий также будут различными, так как они вычисляются по формулам (11), которые имеют разные временные базы для прямого и обратного времен. Здесь имеются в виду статистически значимые отличия, достоверность которых может быть оценена по соответствующему критерию. Аналогичные рассуждения могут быть проведены и относительно ковариационных матриц в целом. Следовательно, в нестационарном приближении в общем случае выполнение равенства (15') не соблюдается даже в одномерном случае. Не выполняется в этом случае и соотношение (15"').
Этот вывод еще в большей степени применим к многомерному случаю, как к объекту более сложной природы. Оценка и учет этого расхождения, по нашему мнению, позволяют смоделировать параметры возможного тренда. Практическая реализация этой идеи заключается в раздельной оценке матриц Db и Dp по не полностью совмещенным временным интервалам. Степень разноса данных зависит как от лага модели l, так и от длины базы исходных данных N.
Итак, после второй итерации получаем матрицы Cb2, Cp2 и вектора остатков B2 и P2, а также оценки их дисперсий Db2 и Dp2 и ковариаций Rbp2, Rpb2. При необходимости выполняем третью итерацию. Для этого рассчитаем по формуле (12) матрицы Cb3, Cp3. При этом вместо векторов Xb2, Xp2 используем вектора остатков Bk2 и Pk2 соответственно, полученные по формулам (14). После применения выражения (13), а затем (14) к этим векторам вновь получаем остатки Bk3 и Pk3, но уже после третьей итерации и т.д. На l-ом шаге получаем матрицы Cbl и Cpl, вектора невязок Bkl, Pkl, оценки их дисперсий и ковариаций.
Совокупность матриц Cbi, Cpi является достаточной для построения оценок фазового вектора, однако практическим неудобством при этом является необходимость использовать не исходные данные, а их разности. Для явного определения параметров моделей (3)-(4) требуется вычислить матрицы Abi, Api по ранее полученным матрицам Cbl, Cpl, т.е. нужно перейти от частных моделей к обобщенным. Для этого воспользуемся рекуррентной формулой Теплица для одномерного случая (Никитин, 1979), которая позволяет по коэффициентам модели (l-1) порядка получить коэффициенты модели l-oro порядка, если известны матрицы Cbl, Cp¡. Обозначим матрицы моделей (l-1) порядка следующим образом: Abi(l-1), Api(l.1),
i = 1,2,___,(/—1). В этом случае матрицы моделей (3-4) l-oro порядка, по аналогии с одномерным случаем,
могут быть вычислены по следующим формулам:
Abl = Cbl — CblAbi(l-1) i = 1,2,> • •Х1-^
Api = Cpl - CplApi(l-i) (членов нет, если l = 1). (16)
После /-ой итерации получаем также ковариационные матрицы ошибок прогноза Dbl и ретроспекции Dp/, достигнутые на /-ом шаге итерации. Совокупность этих данных несет разнообразную и обобщенную информацию относительно синтезируемой модели.
Рассмотрим этот вопрос более подробно.
Будем анализировать эволюцию информации по мере увеличения лага модели. Если принятая к рассмотрению модель адекватна анализируемому процессу, то при увеличении ее порядка / элементы матриц Сы, Cp/ (а также Ab/, Ap/) должны приближаться к нулю, а след матриц Db/, Dp/ - стремиться к постоянным величинам. Большие остаточные значения составляющих матриц Cb, Cp/ свидетельствуют или о малости принятого лага модели, или о нелинейности связей между исходными данными, или о неполном информационном базисе. Последовательно увеличивая лаг модели, можно уменьшить эту неопределенность. Матрицы Db, Dp/ характеризуют немоделируемую часть исходных данных, другими словами, отображают шум модели. Стабилизация оценок этих матриц свидетельствует о том, что достигнута предельная точность выбранной модели, и дальнейшее увеличение ее лага уже не может привести к существенному улучшению оценок, а иногда и ухудшает их. Таким образом, анализ поведения матриц в зависимости от лага модели позволяет определить оптимальное значение порядка разностных уравнений (3-4).
В целом анализ этой информация может быть полезен как при уточнении параметров модели, так и при определении ее порядка. След дисперсионных матриц позволяет судить о достигнутой точности модели как по отдельным координатам фазового вектора, так и о суммарной (усредненной) точности. Внедиагональные элементы характеризуют статистические связи ошибок прогноза фазового вектора и также имеют простую интерпретацию. При увеличении лага модели происходит постепенный учет линейных связей данных при различных временных сдвигах, а невязки становятся все более некоррелируемыми и непредсказуемыми. При адекватной модели в пределе эти невязки должны принимать характер белого шума, т.е. матрицы D должны стремиться к диагональному виду. Числовые значения диагональных элементов этих матриц представляют дисперсию белого шума соответствующей компоненты фазового вектора. Присутствие значимых внедиагональных элементов сигнализирует о наличии нелинейных связей между составляющими временного ряда, или о неполном информационном базисе.
Аналогичным образом можно проинтерпретировать стремление к нулю элементов матриц С, A по мере увеличения лага модели. Известно, что белый шум имеет ковариационную зависимость типа дельта-функции, другими словами, два его различных временных отсчета статистически не связаны, и, следовательно, он не может нести какой-либо прогностической информации, что и отражается в поведении матриц С, A. Рассмотрим эту проблему с несколько другой стороны. Будем считать, что выражения (3-4) представляют модели с постоянными коэффициентами, т.е. предположим, что матрицы Abi, Ap/ не меняются во времени. Это предположение соответствует стационарному приближению, для которого будет выполняться равенство Ab = A'p/. Таким образом, для классического метода наблюдается равноправность, с точки зрения качества прогноза, прямого и обратного времен. Для рассматриваемой методики, как видно из приведенных выше замечаний, эта равноправность должна нарушаться даже для модели с постоянными коэффициентами. Это нарушение будет проявляться в виде некоторого рассогласования аналогичных матриц, построенных для прямого и обратного времен. Следовательно, будут несколько отличаться и соответствующие прогностические модели и качество их работы. Степень рассогласования аналогичных матриц, или качества прогноза, может являться характеристикой нестационарности процесса. Но особенно хорошим индикатором тренда являются нормированные ФАК, имеющие ясный физический смысл и наиболее простые и наглядные для анализа. Поэтому рассмотрим задачу расчета корреляционных матриц по полученной выше информации.
Покажем далее, что знание матриц Ab/, Ap/ позволяет оценить функции автокорреляции фазового вектора. Для одномерного случая Дж. Бургом показано, что эти параметры связаны следующим матричным соотношением (Никитин, 1979):
R(0) R(1) ... R( ) 1 Di
R(1)R(0) ... R( i -1) Ai/ = 0
R(i) R( i-1) .. R(0) A// 0 (17)
где R(i) - значение функции автокорреляции при /-ом сдвиге.
Это выражение позволяет по известной последовательности чисел Ац, I = 1,2,...,/ и дисперсии ошибки прогноза Б/ вычислить значение ФАК в /-ой точке. Обобщим это уравнение на многомерный случай. Для перспективного прогноза аналог выражения (17) будет выглядеть следующим образом
Я(0) Л(-1) ... й(1) Л(0) ...
Я(/) Я(/-1) ...
R(-/) R(-/+1) E Ab/i = Db/ 0
R(0) Ab// 0 (18)
где Е - единичная матрица размера тхт; Я(к) - искомые субматрицы автокорреляций при сдвиге к; 0 -нулевая матрица размера тхт.
Нетрудно показать, что аналогом формулы (17) в случае ретроспективного прогноза является выражение следующего вида
R(0) R(-1) .. R(-D Apu 0
R(1) R(0) . R(-/+1) ap,(/-1),1 = 0
R(/) R(/-1) .. . R(0) E Dpi
(19)
Уравнение (18) позволяет вычислить матрицу R(/), если известны матрицы R(0), R(1), ... , R(/-1) по следующей формуле
R(l) = -TiR(l-i)Ab//, i = 1,2,...,/. (20)
Аналогичным образом уравнение (19) дает возможность определить матрицу R(-l), если известны матрицы R(0), R(-1), ... , R(-/+1) по следующему выражению
R(-l) = -!i R(-i) Ар/.,-/, i = 1,2,.,/. (21)
Как видно из приведенного алгоритма, для начального старта метода требуется знать матрицу R(0), а остальные могут быть получены по формулам (20), (21). Матрица дисперсий R(0) в данном алгоритме вычисляется по стандартным методикам. Практические расчеты корреляционных матриц по алгоритмам (20), (21) совмещены с последовательностью вычисления матриц Ab/, Ар/ по формулам (16), что способствует минимизации затрат общего времени расчетов.
Отметим, что полученные по данному методу ВФАК имеют принципиальные отличия от аналогичных параметров, рассчитанных по классическому варианту. Для стационарного случая в обеих методиках имеет место равенство
Rb(í) = R'p(-T),
Ab(k) = A'p(k). (22)
Рассмотрим теперь нестационарный процесс. Формально и в этом случае для классического метода соотношения (22) будут выполняться, так как соответствующие матрицы получаются простым их усреднением по прямому и обратному временам. В то же время анализ формул (12-14, 16, 18-22) показывает, что для разработанного приближения соотношение (22) не выполняется и, следовательно, будет наблюдаться некоторое рассогласование корреляционных функций, рассчитанных в перспективном и ретроспективном направлениях. Это явление можно объяснить, как несоответствие статистической структуры процесса в прямом и обратном времени. Этот факт кажется почти тривиальным, если принять во внимание, что прогнозируемость информации зависит от целого ряда статистических факторов, которые могут по-разному локализоваться во времени, определяя неравноправность прямого и обратного прогнозов. В частности, одной из причин, во многом определяющих прогнозируемость процесса, является тренд.
С другой стороны, алгоритмы (20), (21) можно рассматривать как прогноз матриц ВФАК на /-ой временной задержке по их значениям на предыдущих сдвигах с помощью параметров линейной модели
АЬ/, Ар/. Так как, по сравнению с классическим вариантом, эти параметры учитывают более тонкие особенности статистической структуры данных, полученные оценки ВФАК являются статистически более обоснованными и точными и, следовательно, обеспечивают более достоверное продолжение корреляционных функций. Известно (Рожков, Трапезников, 1990), что для одномерных векторов и в стационарном приближении этот метод используется для построения спектров. По сравнению с классическими схемами расчета спектра, этот подход обеспечивает большее разрешение по частоте и устойчивость полученных результатов. Учитывая, что спектр и ВФАК имеют взаимно однозначное соответствие, такой результат обеспечивается за счет достоверной экстраполяции корреляционных функций по нескольким начальным и наиболее достоверным ее значениям. Таким образом, предлагаемая методика позволяет строить более достоверные оценки ВФАК с более высоким разрешением и точностью даже при больших лагах модели. Все эти рассуждения могут быть перенесены и на более общий, многомерный стационарный случай.
Мера рассогласования корреляционных матриц для прямого и обратного времен может служить критерием тренда процесса. Вычисление автокорреляций прогноза и ретроспекции по реальным данным показывает, что это рассогласование проявляется в форме систематического смещения этих функций относительно друг друга (даже для одномерного случая). Степень этого рассогласования как раз и связана с систематическим смещением составляющих фазового вектора. Более тщательный анализ этой проблемы показывает, что имеется принципиальная возможность разделить корреляционную функцию на трендовую, "сигнальную" и "шумовую" компоненты, что является шагом к явному определению этих составляющих фазового вектора. Корреляционные матрицы могут быть использованы и для проверки полноты информационного базиса. Здесь под полнотой информационного базиса понимается возможность в рамках линейного подхода по имеющимся в наличии данным построить модель с заданными свойствами по точности. Успешный прогноз может быть построен только для достаточно связанных данных, что может быть проверено по матрицам взаимной корреляции.
В данной статье рассмотрен одношаговый предсказатель, однако метод легко может быть изменен и для прогноза с большей задержкой. Действительно, пусть требуется осуществить прогноз на 5 шагов (5 > 1). Для этого достаточно в качестве стартовой итерации выбрать итерацию с номером /н = 5 + 1 (напомним, что для одношагового прогноза стартовая итерация равна 2). Далее приведенная выше методика может быть использована без изменений.
Отметим также, что стандартные методы решения этой задачи для стационарного случая требуют затрат времени ЭВМ пропорционально (тх/)3 и объемов машинной памяти пропорционально величине (тхI)2. Для матриц Теплица, свойства которых использованы в данном алгоритме, требуются затраты времени пропорционально (тх/)2 и объем памяти пропорционально величине 2хтхI.
4. Заключение
Разработанный метод был протестирован на реальном геофизическом материале. Выполненная проверка показала, что он позволяет моделировать достаточно сложные геофизические процессы. Алгоритмы обладают повышенной устойчивостью даже при порядках моделей (3-4), близких к длине интервала наблюдения. Традиционные методы при этом не работают из-за вырождения обращаемых матриц. Достигаемая точность метода на 15-30 % выше, чем у классических подходов (при одинаковых порядках моделей). Улучшение достигается как за счет учета нестационарности анализируемых процессов, так и применением методически более обоснованных способов анализа информации.
Особенно эффективен алгоритм в сочетании с методом нейроподобного моделирования, где компактность и быстродействие играют решающую роль. Так как при моделировании попутно вычисляются автокорреляционные матрицы, методика позволяет успешно решать многие задачи многомерного спектрального анализа. Более того, как и в одномерном случае (Никитин, 1979), многомерный спектр может быть построен на основе матриц АЬ, Ар.
Важным достоинством метода является принципиальная возможность явного определения трендовой, "сигнальной" и "шумовой" составляющих процесса.
Литература
Клаербоут Дж.Ф. Теоретические основы обработки геофизической информации. М., Недра, 304 е., 1981. Котюк А.Ю., Ольшевский В.В., Цветков Э.И. Методы и аппаратура для анализа характеристик
случайных процессов. М., Энергия, 237 е., 1967. Макузе Ю.И. Алгоритмы для уравнивания геодезических сетей на ЭВМ. М., Недра, с.56-57, 1989.
Никитин A.A. Теоретические основы обработки геофизической информации. М., Недра, 342 е., 1979.
Pao С.Р. Линейные статистические методы и их применение. М., Наука, с.73, 1968.
Рожков В.А., Трапезников Ю.А. Вероятностные модели океанологических процессов.
Л., Гидрометеоиздат, 272 е., 1990. Сильвия М.Т., Робинсон Э.А. Обратная фильтрация геофизических временных рядов при разведке на
нефть и газ. М., Недра, 248 е., 1983. Троян В.М., Соколов Ю.М. Методы аппроксимации геофизических данных на ЭВМ. Л., Ленинградский
университет, 302 е., 1989. Химмерблау Д. Прикладное нелинейное программирование. М., Мир, 534 е., 1975.
