CC BY
82
Моделирование одномерных динамических процессов с целью предварительной обработки результатов
Ю.П. Драница
Судоводительский факультет МГТУ, кафедра высшей математики и программного обеспечения ЭВМ
Аннотация. Разработана теория предварительной обработки одномерных рядов данных, основанная на линейной динамической модели. Показано, что на основе этого метода можно решать широкий круг задач, связанных с предварительным анализом и преобразованием информации, а именно: определение параметров шума и его фильтрация; создание базиса для аппроксимации, сглаживания и экстраполяции функций; надежное оценивание функций автокорреляции; спектральный анализ высокого разрешения и другие проблемы. Теория предназначена для анализа одномерных временных рядов. Возможности методики продемонстрированы на конкретном натурном материале. Разработанная методика имеет хорошее теоретическое обоснование и опирается на физическую сущность изучаемого явления.
Abstract. The theory of the preliminary processing of the one-dimensional sequences of data based on the linear dynamic model has been worked out. It has been shown that the problems of preliminary analyses and information transforming can be solved with the help of this method. The theory is intended for the analysis of one-dimensional time series. The possibilities of the method have been shown on a concrete experimental material. The designed method has a good theoretical substantiation and is based on the physical essence of the studied phenomenon.
1. Введение
При изучении сложных систем и процессов исключительно важное значение имеет предварительная обработка исходных данных, которая включает широкий круг задач, облегчающих последующий анализ и интерпретацию информации. К таким задачам, например, относятся определение шумовых характеристик сигнала и отношения сигнал-шум, фильтрация шумовой компоненты, выбор базисной системы функций, в которой предполагается разлагать данные, достоверная оценка функций автокорреляции (ФАК), экстраполяция сигналов, фильтрация в заданной полосе частот, вычисление спектров с высоким разрешением, сглаживание функций и многие другие вопросы. Причем желательно решать все эти задачи, опираясь на физическую сущность изучаемого явления, а не использовать чисто формальные приемы. Анализ этих проблем показывает, что многие из указанных выше задач могут быть успешно решены в рамках линейного динамического моделирования, которое имеет хорошее теоретическое обоснование и развитый математический аппарат.
2. Общая постановка проблемы
Рассмотренные ниже методики предназначены в основном для анализа дискретно-временных или дискретно-пространственных последовательностей, которые получаются в результате дискретизации сигналов, непрерывно изменяющихся в пространстве или времени. При этом предполагается, что используется равномерная дискретизация. Рассмотрена также задача перехода от дискретного представления моделей к непрерывному.
Рассмотрим уравнение одномерной авторегрессии /-ого порядка
Хк = ХАХм+Ъ i = 1,2,...,/, (1)
где Хк - одномерный ряд; к = 1,2,...; Ai - коэффициенты авторегрессии; / - лаг модели (максимальное число удерживаемых в уравнении (1) членов временного ряда для прогноза его k-ого элемента); sk - ошибка прогноза.
Поставим в соответствие модели (1) разностное уравнение с постоянными коэффициентами (Корн, Корн, 1977)
EAXk-i = fk), i = 0,1,., /, Ao= 1, (2)
где f(k) - внешние возмущения процесса. Рассмотрим однородную часть уравнения
ZAXk-, = 0. (3)
Как известно (Корн, Корн, 1977), решение однородного уравнения (3) имеет следующий вид:
X = С1(Я1)к+ С^+.+ОД^, (4)
где С, - некоторые константы; / = 1,2,...,/; Я2, ...Д - корни характеристического уравнения
ЛоЛ' + Л^'"1 + ...+ Л\ = 0, (5)
если только все его / корней различны. Если некоторый корень, например 1, имеет кратность т, то соответствующий член в решении (4) будет иметь вид:
(С^^+.-.+г-СтКАА (6)
При действительных коэффициентах уравнения (5) двум простым комплексно-сопряженным корням Я = ехр(о±/®), можно поставить в соответствие решение
ехр(ак)[Лсо8(к®)+В8ш(к®)], (7)
где А, В - неизвестные константы; к - дискретный параметр времени.
Таким образом, любому уравнению авторегрессии (1) может быть поставлено в соответствие разностное уравнение (3) с решениями (4-7). Известно, что решение (4) будет устойчивым, если выполняется неравенство А< 1, а решение (7) - при а< 0. В этом случае решение будет представлять сумму затухающих экспонент, или затухающих по экспоненте косинусоид различных частот с коэффициентами затухания а.
Будем считать, что уравнения (1-2) являются результатом дискретизации соответствующих непрерывных дифференциальных уравнений по времени или другой упорядоченной переменной, т.е. принимаем, что к есть дискретный параметр времени. Упорядочим решения этих уравнений по возрастанию частот и примем их в качестве системы базисных функций:
Yl(k) = С1ехр(«1 к),., Yj = С]ехр(а]к), Yj+l(к) = exp(aj+lк)[Лj+lcos(aj+lк)+Bj+lsm(щ+lк)],..., (8)
^(к) = exp(a\к)[ЛlCOs(&lк)+Bj+lSm.(a)lк)],
где нижний индекс обозначает номер функции в базисе, а аргументом является параметр дискретного времени; А, В, С - неизвестные коэффициенты; 3 - число базисных функций с нулевыми частотами.
Сформированный базис будем использовать в дальнейшем для различных преобразований информации. Отметим, что для этих целей могут быть использованы достаточно разнообразные системы функций, например, полиномиальные, тригонометрические, полиномы Чебышева и другие. Однако в большинстве случаев эти базисные системы имеют чисто формальный характер, мало связанный с существом обрабатываемой информации. Например, в разложении Фурье аппроксимирующие частоты обусловлены периодом наблюдений, которые никак не связаны с физической сущностью отображаемых процессов, полиномиальная аппроксимация также в большинстве случаев не имеет какого-либо физического содержания. Напротив, данная базисная система (8) имеет четкое физическое обоснование, обусловленное представлением процесса линейной динамической моделью (1) с решениями (4-7).
Анализ результатов аппроксимации реальных геофизических данных показывает, что это разложение имеет ряд особенностей. Отметим некоторые из них. Прежде всего наблюденные частоты не кратны как периоду наблюдений, так и между собой, хотя при этом может наблюдаться несколько компонент, имеющих нулевые частоты, но разные коэффициенты затухания. В целом экспоненциальные составляющие имеют достаточно высокие коэффициенты затухания и, следовательно, могут иметь значимые амплитуды только вблизи начала координат. Коэффициенты затухания для косинусных составляющих значительно меньше. Как правило, частоты располагаются на числовой оси неравномерно, причем с возрастанием частоты коэффициенты затухания в основном убывают.
Сформированный таким образом базис позволяет достаточно эффективно решать ряд прикладных задач, таких, например, как экстраполяция функций, оценка шумовой компоненты ряда, его разложение по системе функций и другие. Рассмотрим сначала задачу разложения ФАК в этом базисе.
3. Алгоритмы анализа информации. 3.1. Предварительная обработка ФАК
Пусть требуется разложить ФАК временного ряда Х(к) по данному базису. В работе (Драница, 2000) излагается корректный алгоритм расчета функции автокорреляции и уравнения регрессии (1). Воспользуемся этим алгоритмом и рассчитаем ФАК временного ряда с заданным лагом Т. Полученную таким образом функцию обозначим как Я(к). Сформируем матрицу значений базисных функций,
принимая единичными значения амплитуд в базисе (8), т.е. полагаем C1=Ai=Bl=\. При предположениях указанная матрица будет иметь следующий вид:
' 71(0) 72(0) ... Ym(0) ' ^(1) 72(1) ... 7т(1)
A =
Yi(T-l) 72(Г-1) ... YJT-1) Yi(T) Y2(T) ... Ym(T)
где m - общее число функций, которые будут использоваться при разложении; T- максимальная задержка ФАК. Разложение ФАК по этому базису в матричных обозначениях будет иметь следующий вид:
R = AZ, (10)
где Z - вектор коэффициентов разложения (весовые коэффициенты); R - вектор коэффициентов ФАК, состоящий из Тэлементов.
Сумму квадратов S невязки между разложением (10) и анализируемой ФАК можно представить следующим образом:
S(Z) = (R-AZ)'(R-AZ), (11)
где значком " ' " обозначена операция транспонирования.
Проблема состоит в определении оптимальных (по некоторому критерию) значений вектора Z. При m < T аппроксимацию (10) можно рассматривать как типичную задачу метода наименьших квадратов (МНК) с критерием оптимальности в виде минимума суммы квадрата невязки (11). Нам необходимо найти такие оценки вектора Z параметров модели (10), которые удовлетворяли бы этому условию. Для раскрытия экстремума функций S(Z) воспользуемся правилами векторного дифференцирования (Pao, 1968), после использования которых находим необходимое условие экстремума
dS/dZ = - R'R+R'AZ=0,
откуда имеем
AAZ = R'A. (12)
Условие минимума квадратичной формы S означает, что матрица
d2S/dZ2 = A'A
должна быть положительно определена.
Выражение (12) представляет собой матричную запись так называемой нормальной системы линейных уравнений. Если матрица A'A неособенная, то система (12) имеет единственное решение следующего вида:
Z = (A'A)-1R'A, (13)
где символом " -1" обозначена операция вычисления обратной матрицы.
Полученное решение (13) позволяет представить ФАК следующим образом:
R(k) = E,exp(aIk)(AIcos(®IkI)+BIsin(®Ik)), i = 1,2,.,m, k = 1,2,.,T0. (14)
Соотношение (14) можно представить в несколько ином виде:
R(k) = £ iC iCOs(®k+^¿), (15)
где % = arctg(Bi/Ai) - начальная фаза компоненты разложения; С i = (A^+B,2)0'5 - обобщенная амплитуда.
Отметим, что параметр T0 может удовлетворять как условию T0 < T, так и T0 > T. В первом случае производится интерполяция сигнала, а во втором - его экстраполяция, или прогноз. Из выражений (14), (15) следует, что Z является вектором весовых коэффициентов, с которыми базисные функции входят в разложение ФАК. Для нулевых частот (щ = 0) коэффициенты Bi в формуле (14) равны нулю и, следовательно, начальные фазы в разложении (15) также будут нулевыми (pi = 0). Следовательно, нулевые частоты являются частным случаем косинусного представления и входят в разложение (14-15). Очевидно, что данная аппроксимация может быть использована и для любой другой функции ф(к), имеющей ровно T отсчетов. Для этого достаточно подставить эту функцию в уравнение (13) вместо вектора R и получить решения (14-15).
Отметим, что матрица А необязательно должна включать все базисные функции, а может представлять некоторое (весьма произвольное) их подмножество. Следовательно, разложение можно выполнять как по полному, так и по некоторому частичному базисам. Обозначим через In, I1, I2,
соответственно, полный и какие-либо два частных базиса, а через Sn, Sb S2 - разложение ФАК по ним. Введем третий частный базис I3 = Ii+I2 и обозначим разложение ФАК по нему S3. Тогда, если базисы I и I2 не имеют общих функций, то, вследствие линейности операций разложения, будем иметь следующее равенство
S3 = Si+ S2. (16)
Соотношение (16) выражает принцип суперпозиции решений, следующий из теории линейных систем. Рассмотрим точность аппроксимации (14). Из самых общих соображений ясно, что она зависит, в первую очередь, от точности решения уравнения (13), числа удерживаемых в разложении базисных функций и сложности разлагаемой функции. Экспериментами установлено, что для функции типа ФАК точность решения уравнения (13) определяется выбранным численным методом решения системы уравнений (13), а также разрядностью используемой арифметики ЭВМ и слабо зависит от сложности аппроксимируемой ФАК. В среднем относительная точность решения уравнения (13) составляет более 99.999 %, т.е. в практическом плане этими ошибками можно пренебречь. Этот интересный факт вытекает из того, что исчерпывающей информацией о модели (1) является ковариационная функция временного ряда X(t). В частности, этой информации достаточно для определения коэффициентов регрессии (1), т.к. существует прямая связь между ними и функцией ковариации (Губанов, 1997).
Число членов разложения, которое требуется использовать для достижения требуемой точности восстановления аппроксимируемой функции, существенно зависит от ее гладкости. При переходе от преимущественно низкочастотных функций к более высокочастотным в общем случае требуется использовать большее число членов разложения (14-15). Экспериментами установлено, что для восстановления ФАК с точностью 95-99 %, в расчет необходимо принимать 3-7 (в зависимости от ее гладкости) базисных функций при лаге T = 50. Количественно точность аппроксимации может быть определена следующим образом. Обозначим через Zm, Si - сумму квадратов коэффициентов разложения, соответственно, по полному и какому-либо частному базисам. Тогда, вследствие принципа суперпозиции, точность аппроксимации ФАК может быть представлена следующим выражением:
4 = h/Lm . (17)
Таким образом, имея разложение ФАК по полному базису, по соотношению (17) можно затем осуществить ее аппроксимацию с заданной точностью. Очевидно, что приближение функции с наперед заданной точностью может быть осуществлено по нескольким базисам, т.е. в общем случае не является однозначной операцией. Следовательно, для однозначного приближения функции необходима дополнительная информация.
Отметим, что важным результатом представленного выше алгоритма аппроксимации ФАК является то, что функция приближается достаточно простыми аналитическими выражениями, которые облегчают решение многих задач, например, расчет спектра, физическую интерпретацию результатов и другие. Кроме того, аналитическое представление позволяет очень просто осуществить переход к непрерывному времени (т.е. перейти от конечных разностей к обыкновенным дифференциальным уравнениям). Известно, что теория обыкновенных дифференциальных уравнений тесно связана со многими приложениями, такими, например, как линейные динамические системы, электро- и радиотехника и т.п. Выводы этой теории имеют строгую и простую физическую интерпретацию и позволяют наиболее адекватно описать и интерпретировать многие процессы и системы. Интерпретация результатов развиваемого метода в рамках этой теории излагается ниже.
3.2. Оценка параметров шума
Рассмотрим возможности метода для оценки параметров шума сигнала, описываемого уравнениями (1) или (2). Обозначим спектральную плотность сигнала функцией S(®). Как известно, спектральная плотность белого шума не зависит от частоты, т.е. S(®) = const и имеет ФАК типа дельта-функции на нулевом сдвиге. На реальной ФАК эта особенность спектра проявляется в виде скачка (или излома) функции в начале координат. Величина этого скачка равна дисперсии некоррелированного белого шума. Анализ структуры базисных функций показывает, что такому поведению ФАК соответствуют компоненты с нулевой частотой (чисто экспоненциальные члены) и большим коэффициентом затухания. Такие базисные функции могут вносить наибольший вклад в разложение только на сдвигах ФАК, близких к началу координат (1-3 отсчета), далее эта составляющая становится пренебрежимо малой по сравнению с компонентами с ненулевыми частотами. Такое поведение компонент разложения с со= 0 можно проинтерпретировать как проявление слабо коррелированного белого шума. Примем эти компоненты аппроксимации функции за ФАК слабо коррелированного белого шума на нулевых частотах.
Из априорных соображений обычно бывают известными наибольшие частоты /р), которые для данного процесса имеют определенный физический смысл, а все частоты на ФАК с / > /р следует рассматривать, как проявление коррелированного (цветного) шума. Обозначим числом М1 - номер наибольшей базисной функции, для которой еще выполняется условии/</тр. Выполним теперь алгоритм (8-10) и разложим ФАК в базисе с номерами функций в диапазоне 3+1^Ы1. В результате получим разложение функции, очищенное от слабо коррелированного белого шума. Разность между исходной функцией и ее аппроксимацией даст ФАК белого шума. Иногда бывает необходимым выделить одиночную компоненту сигнала с частотой, лежащей в интервале/с. (0/р], например помеху, связанную с частотой питающей сети. Для этого достаточно исключить функцию с требуемой частотой из базиса. Другим способом выделения шума является формирование частного базиса, компоненты которого связаны с шумом, и последующим разложением ФАК по этой системе функций.
Таким образом, аппроксимация ФАК процесса и частотный анализ этого разложения позволяет оценить корреляционные (а следовательно и спектральные) характеристики действующих шумов и выделить их из коррелограммы. Эти параметры шума объективны в том смысле, что они следуют из решения уравнений (1) или (2) и являются отражением физической сущности протекающих в системе процессов. Важным достоинством этого метода является то, что идентифицируемый при этом шум может являться цветным. Решение этой задачи классическими методами встречает большие затруднения и нередко носит чисто формальный характер.
3.3. Задача надежной оценки, сглаживания и экстраполяции ФАК
При анализе коротких рядов данных часто возникает задача надежной оценки корреляционных функций. Известно, что достоверная оценка ФАК возможна только для сдвигов, в 5-10 раз меньших, чем длина анализируемого ряда (Котюк и др., 1967). При больших сдвигах оценка становится неустойчивой и может не удовлетворять нужным свойствам. Одним из важнейших требований, которому должна удовлетворять эмпирическая оценка ФАК, является ее положительная определенность. Только в этом случае она позволяет адекватно рассчитывать спектр и может быть использована в МНК. Известно, что если ФАК не является положительно определенной, то спектр процесса может иметь отрицательные компоненты, а оценки параметров линейной модели не будут обладать оптимальными свойствами. Обычно эта проблема решается следующим образом. Исходя из длины анализируемого ряда, находят несколько достоверных оценок ФАК, которые затем экстраполируются на требуемый период времени. Экстраполяция осуществляется с соблюдением требований положительной определенности функции.
Аналогичный подход применен и в разрабатываемом методе. Используя методику (Драница, 2000), оценим несколько начальных и наиболее достоверных отсчетов ФАК, а также построим систему базисных функций по методике, описанной выше. Разложим достоверную часть ФАК по данному базису (т.е. получим веса базисных функций Л). Далее, имея вектор весов Л, легко получить экстраполяцию ФАК на требуемый лаг по формулам (14, 15). Существенным достоинством данного подхода является то, что здесь применено аналитическое продолжение функции, что может оказаться весьма важным при дальнейшем использовании ФАК. Отметим, что такая процедура автоматически гарантирует получение положительно определенной ФАК. Это следует из того, что спектр суммы затухающих косинусоид всегда положителен. В то же время ФАК при небольших сдвигах имеет, как правило, хорошую гладкость, следовательно, эта процедура позволяет осуществить также некоторое сглаживание этой функции.
Другой способ сглаживания ФАК заключается в ее разложении по какому-либо частному базису. Для этого необходимо соответствующим образом сформировать матрицу А (выражение (9)) и получить аппроксимацию ФАК по этому частному базису. Вполне очевидно, что этот базис должен содержать наиболее низкочастотные базисные функции (за исключением компонент с со, = 0). Сглаживание может сопровождаться также экстраполяцией функции. Очевидно, что этот способ аппроксимации (при соответствующем выборе частного базиса) можно рассматривать как своеобразный способ полосовой фильтрации ФАК.
3.4. Разложение временного ряда по системе базисных функций
Рассмотрим теперь задачу разложения временного ряда по базису (8). Предположим, что имеется одномерный ряд Х(к), для которого уже построена система базисных функций, оценена ФАК (Я(к)) и его шумовой компоненты (Яш(к)). Будем производить разложение временного ряда в построенном базисе. Обозначим составляющие будущего разложения следующим образом:
Х1(к)Х£к),...Жк), 0<к<Тэ, (18)
где I - число векторов, входящих в базис; нижний индекс обозначает номер компоненты разложения; Т0 - длина временного ряда.
Будем считать, что вектор данных Xизвестен, и требуется оценить вектораХ1, Х2,..., Х1. Введем обозначение Б = X, и предположим, что оценки этих сигналов могут быть представлены в классе линейных оценок. Пометим оценку вектора сигналов 8 значком " Л " и примем
НХ, (19)
где Н - некоторая матрица размером I х Т0.
Иначе говоря, каждая компонента вектора X, аппроксимируется линейной комбинацией вектора X. С другой стороны, исходную функцию Х(к) можно интерпретировать как суперпозицию пока неизвестных сигналов ХьХ2,...Хь Таким образом, мы полагаем, что исходный ряд Х(к) представляет некоторое множество "измерений", а компоненты (18) - некоторые "сигналы", которые необходимо выделить из этой смеси. Предполагается, что все эти векторные величины центрированы, т.е.
Е{Х}= 0; Е{Б} = 0, (20)
где Е{} - обозначена операция вычисления математического ожидания. Тогда выражение е = 5^-Х будет представлять вектор ошибок оценивания, со следующей ковариационной матрицей:
2ее = Е{вв'} = Е{(5Т-Х)(5Г-Х)'}. (21)
Диагональные элементы этой матрицы представляют собой дисперсии компонент оцениваемого вектора Б.
о2= ЕК^-Х)2}. (22)
Согласно общей теории статистического оценивания, наилучшей оценкой сигнала Б по данным Х является линейная оценка вида (19) с минимальной дисперсией. Усредняя по реализациям выражение (19) и принимая во внимание (20), имеем
Е{БЛ} = НЕ{Х} = 0 = Е{Х}.
Поэтому оценка (19) является несмещенной при любой матрице Н.
Определим эту матрицу таким образом, чтобы дисперсии ошибок оценивания (21) были минимальны. Для этого распишем выражение (21) более подробно
Qee = Е{(НХ-Б)(НХ-Б)'} = НЕ{(ХХ}Н'-Е{БХ}Н-НЕ{ХБ'}+Е{ББ'}= Н2ххН'-&хН'-Н2Х8+&8, (23)
где Цхх, Ццц - автоковариационные матрицы векторов Х и Б соответственно; Qsx - матрица взаимных ковариаций между Б и Х.
Предполагается, что эти матрицы имеют полный ранг (говорят, что кхш матрица А имеет полный ранг, если гапк(А) равен меньшему из чисел к и т). В этом предположении матрица Охх имеет обратную матрицу, и тогда непосредственной проверкой можно убедиться, что выражение (23) эквивалентно следующему:
2ее= 288-28х(2хх)-12х8+[Н-28х(2хх)-1]2хх[Н-28х(2хх)-1]'. (24)
Обозначим А = Qss-QSx(Qxx)'lQxS и В = [Н-28х(2хх)-1]2хх[Н-28х(2хх)-1]'. Тогда выражение (24) будет иметь следующий вид:
2ее = А+В.
Матрица А не зависит от матрицы Н и поэтому одинакова при любых оценках (19). В работе (Губанов, 1997) показано, что значения диагональных элементов матрицы В положительны, следовательно, минимальное значение предыдущего выражения и ошибки (24) наблюдаются при
Н = &х(2хх)-1 (25)
и выражение (18) будет выглядеть следующим образом
8Л= &х(2хх)-1Х. (26)
Конкретизируем выражения (25-26) в контексте решаемой проблемы. Обозначим компоненты разложения ФАК Я(к) временного ряда X как Ль Л2,..., Л/, и определим матрицу Qxx и вектор Qsx следующим образом:
-Л(0) Л(1) ... Я(Т0) Я(1) Я(0) ... Л(Т0-1)
Qxx =
Л(Т)-1) Л(Т0-2) Я(Т0) Л(Т0-1)
Qx =
Я,(0) Я,(1)
Я,(1) Я,(0)
ВДг1) Л,(Т0-2) Я^) Л,(Т0-1)
Я(1) Я(0)
Яг(Т0) я,(Т0-1)
Я,(1) Я ,(0)
(28)
Ранее было сделано предположение, что матрица Qxx имеет полный ранг. Однако это требование не является жестким, так как в противном случае весь необходимый аппарат можно выполнить на основе сингулярного приближения (Троян, Соколов, 1989).
Рассмотрим выражение (26) более подробно. Множитель QsX(QXx)~l представляет матрицу размера Т0 х Т0, элементы которой представляют некоторый цифровой фильтр. Следовательно, преобразование (19) можно рассматривать как некоторую фильтрацию данных. В результате этой фильтрации будет выделен сигнал с заданными частотными свойствами (с частотой ®и затуханием а). В частности, соотношение (19) пригодно для вычленения шумовой компоненты сигнала, для этого необходимо вектор представить через ковариацию шума Яш.
Самой характерной особенностью этого фильтра является то, что он способен выделять составляющие сигнала, имеющие характер затухающих экспонент и затухающих косинусоид с заданными корреляционными параметрами, не нарушая остальных составляющих сигнала. Другими словами, фильтр способен селектировать компоненты решения (4-7). Таким образом, полученные разложения сигнала имеют определенный физический смысл. Заметим, что при выводе этих соотношений не делалось никаких предположений о числе "сигналов" и количестве данных, поэтому они справедливы при любой их величине. Это означает, что мы можем с помощью формулы (26) оценить любое (в принципе неограниченное) количество сигналов Хь Х2, ... по конечному числу данных, если, конечно, для них известна ковариационная матрица 0«. В частности, можно выделить некоторые комбинации решения. В контексте решаемой задачи разложение стало возможным вследствие факторизации ФАК по системе базисных функций и принципа суперпозиции. И в этом случае появляется возможность оценки вектора 0«, компоненты которого представляют, по сути дела, корреляцию, рассчитанную по каждой базисной функции.
3.5. Фильтрация в заданной полосе частот и переход к непрерывному времени
Рассмотренная выше методика фильтрации позволяет, по сути дела, выделять из исходного сигнала всего одно решение вида (4-7) с определенными частотными свойствами. Рассмотрим возможность выделения из смеси группы таких решений, т.е. осуществим фильтрацию в некоторой полосе частот. Пусть необходимо выделить составляющие сигнала в полосе частот а с [<% < ю< ®,2]. Сформируем частный базис /ч таким образом, чтобы все входящие в него базисные функции удовлетворяли заданному выше диапазону частот. Введем следующее обозначение:
= Е,Е{ХХ) = ЕЛ, г = /1,/ 1+1,.,,2.
При данном обозначении выражение (26) преобразуется к следующему виду:
Хч-а = е^й«)-1*. (29)
Очевидно, что сигнал в левой части уравнения (29) представляет суперпозицию сигналов с частотами в заданном диапазоне. Это следует из линейности используемых операций и принципа суперпозиции. Аналогичным образом может выделяться сигнал, представляющий более сложную комбинацию частот.
Рассмотрим непрерывный аналог уравнения (1), т.е. дифференциальное уравнение Т-ого порядка. Тогда процесс, который описывается этим уравнением, можно представить в виде нескольких независимых подсистем, каждая из которых возмущается белым шумом ¿(/) и описывается следующим стохастическим дифференциальным 2-ого порядка с постоянными коэффициентами (Корн, Корн, 1977)
а0¿^х/й^+а^4х!Ж+а2х = Ь0е(/),, = 1,2,.,/.
(30)
Число уравнений в (30) равно количеству значимых компонент разложения ФАК процесса (/). Выходом такой системы является суперпозиция выходов каждой подсистемы. Известно, что решения уравнения (30) имеют корреляционные функции вида (15). Такая связь позволяет, в частности, восстановить коэффициенты каждого уравнения системы (30). Запишем уравнение (15) в несколько иной форме
Щ) = £г(о"/со8^)ехр(-а;0со8(й>/+^), г = 1,2, ...,/, (31)
где (Гг - дисперсия белого шума составляющих случайного процесса Хг.
В этих обозначениях постоянные коэффициенты системы (30) будут выражаться следующим
образом
аг = аи/й0г, а>г = [4а'0г«,2г-(а1г)2]а5/(2а,0г), Ь0\ = (^ш($-^)|)а5, фг = а1ап(агМ). (32)
Обычно полагают, что а0г = 1, ив этом случае система (32), с учетом соотношения (31), имеет два уравнения и два неизвестных и, следовательно, позволяет однозначно определить коэффициенты уравнения (30). Более того, рассмотренные выше решения (4), (6) также являются решениями однородного уравнения (30), которые называются собственными колебаниями.
Таким образом, разработанная методика позволяет не только иметь физически обоснованное разложение процесса, но и получить и идентифицировать математические модели анализируемой системы для непрерывного времени.
Аналитическое представление ФАК процесса дает возможность легко рассчитать и его спектр.
3.6. Гармонический анализ
Рассмотрим другой способ аппроксимации данных. Иногда требуется разложить функцию по системе гармонических функций, имеющих некратные частоты, т.е. получить разложение вида
^ХаДю) = а0 + Еагсо8(®/)+Ьг5ш(®/), г = 1,2,.,/, (33)
где а, Ь - вектора коэффициентов разложения; © - вектор действующих частот.
Аппроксимация (33) довольно часто используется в приложениях, и ее решение относится к проблеме выявления скрытых периодичностей. Другими словами, она относится к задаче распознавания спектральной структуры реальных процессов по результатам их измерений. Соотношение (33), в частности, можно рассматривать как моделирование тренда процесса, основанного на полигармоническом представлении (Высоцкий и др., 1975). Задача выявления скрытых периодичностей будет полностью решена, если вычислены параметры аг, Ь г и т. При этом частоты необязательно соизмеримы или кратны друг другу. Этим, в частности, данная постановка отличается от разложения Фурье, предполагающего факторизацию функции по системе кратных частот с основной гармоникой, обратной периоду наблюдений.
Если действующие в процессе частоты априорно известны, то задача (33) преобразуется в типичную линейную, решение которой вполне тривиально. В противном случае эта задача становится существенно нелинейной. Как и ранее, будем определять действующие в процессе частоты с помощью разностного уравнения следующего вида
Хк-1 + ШХк-/+г + АХк-/-) = 0, г = 1,2,.,/. (34)
Таким образом, уравнение (34) представляет симметричное разностное уравнение. Регрессионным аналогом этого соотношения является одношаговый предсказатель, использующий для оценки значений вектора в (к-1) точке / прошлых и будущих значений функции. Следовательно, это уравнение является физически нереализуемым сглаживающим фильтром. Характеристический полином выражения (34) имеет следующее представление
А1Л21+.+А1Л1+1+А,1+А1Л1-1+.+А1 = 0. (35)
Как видно из (35), полученный полином обладает симметричными коэффициентами, причем А0= 1. В большинстве случаев корни этого полинома равны 1. И если к тому же эти корни являются комплексно-сопряженными, то решение будет представлять сумму гармонических колебаний
7(к) = А ¿со5(к®,)+Вг«т(к®,), (36)
где Аг, Вг - неизвестные константы.
В случае, если корни не являются комплексно-сопряженными, решение будет представлять некоторую константу, которая в разложении (33) учитывается свободным членом. После определения действующих в процессе частот, как и в предыдущем случае, может быть сформулирована система базисных функций, представляющая набор синусоидальных функций идентифицированных частот. Следовательно, может быть проделана вся обработка материала, предусмотренная для случая решений в виде суммы затухающих по экспоненте синусоид.
4. Апробация методики на натурных данных.
4.1. Некоторые особенности принятой схемы факторизации данных
Основной проблемой при использовании разработанных выше способов факторизации натурной информации является задача определения лага аппроксимирующей регрессионной модели. Рассмотрим эту проблему сначала для гармонического разложения функций. В этом случае лаг модели должен соответствовать числу реально имеющихся в данных гармонических компонент. В идеальном случае (при отсутствии помех и вычислительных ошибок) при этом должна наблюдаться нулевая невязка между исходными данными и их аппроксимацией. Поэтому для данного случая алгоритм определения порядка модели весьма прост: постепенно увеличиваем лаг модели и рассчитываем невязку между натурными данными и их аппроксимацией по модели с принятым лагом. Момент достижения нулевой невязки и принимается за требуемый лаг аппроксимирующей модели. И в этом случае решение характеристического уравнения будет в точности соответствовать действующим в системе частотам.
В реальных условиях в натурных данных всегда присутствует некоторый шум наблюдений, а все вычисления выполняются с ограниченной точностью за счет конечной точности арифметических операций в ЭВМ. Наличие в данных шума приводит к появлению дополнительных корней решения, не связанных с имеющимися физически обоснованными решениями. Кроме того, выборка данных всегда ограничена, что приводит к погрешностям в определении параметров регрессии. Более того, нет никакой уверенности, что используемая модель является вполне адекватной натурным данным (если для этого нет реальных физических предпосылок). В этих условиях речь может идти только о некотором более или менее удачном решении данной задачи. Поскольку наилучший порядок модели заранее неизвестен, на практике обычно приходится испытывать несколько ее порядков. Полученные таким образом данные подвергают проверке тем или иным критерием с целью окончательного выбора порядка модели.
При аппроксимации данных по системе затухающих по экспоненте косинусоид ситуация несколько иная. Даже если исходные данные неадекватны принятой линейной модели с постоянными коэффициентами, увеличивая порядок модели, можно с какой угодно точностью аппроксимировать информацию в рамках этого подхода. Однако проблема существенно усложняется в условиях ограниченной выборки данных и действия шумов. В этом случае слишком высокий порядок модели приводит к большим ошибкам в оценивании ее параметров, слишком малый - к неточному описанию процесса. Следовательно, для решения этой проблемы по ограниченным и зашумленным данным требуется некоторый компромисс между точностью и адекватностью модели. В большинстве случаев эта задача решается или на основе специальных критериев, или чисто эмпирическим путем. К наиболее распространенным критериям относятся две целевые функции Акаике (Губанов, 1997).
Первый критерий основан на вычислении окончательной ошибки предсказания (ООП). Обозначим ошибку линейного предсказания линейной модели Р-ого порядка символом рр. Тогда первый критерий Акаике определяется выражением
ООП(£>)= ^[^+(р-1))/(^(р-1))], (37)
где N - число отсчетов данных.
Заметим, что член в квадратных скобках растет с увеличением порядка р, характеризуя тем самым неопределенность в оценке ошибки предсказания, обусловленную конечностью выборки данных. С другой стороны, при увеличении лага модели наблюдается уменьшение ошибки предсказания Рр. Выбирается такое значение порядка, при котором выражение (37) минимально. Второй критерий Акаике получил название информационного (ИКА), так как он основан на минимизации некоторой теоретико-информационной функции. В предположении нормальности статистик процесса, ИКА критерий определяется выражением
ИКА(^) = МиОсд + 2р. (38)
Член 2р в выражении (38) характеризует плату за использование дополнительных коэффициентов авторегрессии. Выбирается порядок модели, доставляющий минимум функции (38). При N ^ да эти два критерия асимптотически эквивалентны. Для идеальных авторегрессионных процессов эти критерии дают отличные результаты, однако применительно к реальным данным эти критерии дают заниженные оценки порядка модели.
4.2. Результаты моделирования по натурным данным. 4.2.1. Анализ ФАК процесса
Иллюстрацию применения разработанной теории проведем на временном ряде, представленном на рис. 1 (кривая 1). Эти данные являются отсчетами некоторого геофизического параметра, снятыми через равные интервалы времени (всего проведено N = 201 измерение). Для простоты будем считать, что дискретность отсчетов (АТ) равна единице. В этом случае частоты ю выражаются в безразмерных
единицах, а истинная частота в герцах / выражается следующим образом / = 1/(л®ЖДТ). Будем также считать, что ФАК этого ряда с достаточной точностью может быть оценена при лаге Т = 0.15Ж, т.е. принимаем Т = 30. Согласно работе (Драница, 2000), построим авторегрессионную модель этого ряда и рассчитаем его ФАК (рис. 2, кривая 1). Анализ исходных данных (рис. 1) показывает, что процесс имеет хорошо выдержанную периодическую компоненту, что отчетливо проявляется и на ФАК процесса. Анализируемая функция имеет нулевое среднее и возрастающую со временем дисперсию. Эта нестационарность порождает неэквивалентные регрессии в прямом и обратном временах. Следовательно, будут отличаться и параметры решений, полученные по этим уравнениям. В проведенном ниже анализе использованы частоты и коэффициенты затухания, рассчитанные по обоим уравнениям регрессии.
Решение уравнения (5) имеет 15 корней (гпк = 15): два корня дают решение в виде затухающих экспонент, а остальные являются комплексно-сопряженными. Идентифицированные частоты (только положительные) и коэффициенты затухания системы базисных функций представлены в табл. 1.
Рис. 1. Разложение данных на сигнальную и шумовую составляющие: исходный ряд (1); сигнальная часть данных для базиса 3-8 (2); шум наблюдений (3)
Рис. 2. Функция авторегрессии процесса (1) и её приближение по различным базисам 3(2), 3-4(3), 3-5(4)
Рис. 3. Относительные веса гармоник в разложении ФАК
Частоты в табл. 1 даны в долях радиана, а коэффициенты затухания - в абсолютных единицах. В процессе наблюдается две затухающие экспоненты и тринадцать затухающих по экспоненте синусоид различных частот. Отметим, что максимальная частота, которая еще может присутствовать в процессе, равна п (частота Найквиста).
Аппроксимация функции по полному базису (табл. 1) осуществляется с очень высокой точностью (коэффициент корреляции между ФАК и ее представлением составляет R = 0.9999). Обозначим через квадрат i-ого коэффициента разложения. Тогда, согласно выражению (17), относительный вес гармоники в разложении будет определяться следующим образом:
Лг = Ъ/Lm, i = 1,2,.. .,m*. (3 9)
Относительный вклад гармоник разложения в зависимости от их порядкового номера изображен на рис. 3. Анализ этого рисунка показывает, что наибольший вес в разложении имеют около 3-5 базисных функций. Причем наиболее информативными являются низкочастотные компоненты разложения. Относительные веса некоторых частных базисов, построенных по соотношению (17), представлены в табл. 2.
На рис. 2 представлены аппроксимации ФАК, осуществленные по различным наборам базисных функций табл. 2. Как видно из этих изображений, наибольшее расхождение аппроксимирующих функций наблюдается на малых сдвигах. На больших сдвигах ФАК, построенные по разным базисам, отличаются значительно меньше. Для нулевого сдвига наблюдается излом ФАК, обусловленный наличием в системе помехи, имеющей характер белого шума. Этот шум представляет, с одной стороны,
две затухающие экспоненты, а с другой - наиболее высокочастотные составляющие процесса (табл. 1). Исходя из анализа этой информации можно предположить, что информационная часть данных заключена в частотном окне 0.3071 < а < 0.9186, т.е. принимаем остальные частоты за слабо коррелируемый шум. Изображения информационной части ФАК и ФАК помехи представлены на рис. 4.
Таблица 1. Параметры системы базисных функций
Частоты (рад.) 0.0 0.0 0.3071 0.4026 0.5964 0.7470 0.9186 1.116
Коэф-ты затухания 0.09295 13.35 0.07027 0.04805 0.03329 0.06169 0.07589 0.06161
Частоты (рад.) 1.332 1.485 1.679 1.873 2.344 2.634 2.861
Коэф-ты затухания 0.06621 0.08362 0.06879 0.05252 0.09890 0.05841 0.1007
Таблица 2. Веса некоторых частных базисов
Номера гармоник базиса 1 - 2 3 3 - 4 3 - 5 3 - 6 3 - 7 3 - 8 9 - 15
Относительный вес 0.007 0.622 0.817 0.927 0.952 0.977 0.981 0.012
Рис. 4. ФАК процесса (1), её информационная часть (2) и ФАК шума (3)
Рис. 5. ФАК при сдвиге Т0=30 (*--*), экстраполяция ФАК на Т0=60 (—) и непосредственно рассчитанная ФАК с лагом Т0=60 (_)
Рис. 6. Разложение сигнала по системам базисных функций: исходный ряд (1); базис 3(2); базис 3-4(3); базис 3-5(4)
этой информации показывает, что фак шума не представляет собой дельта-функцию на нулевом сдвиге, т.к. имеет значимые отсчеты на больших запаздываниях. Следовательно, наблюденный шум процесса следует отнести к категории коррелируемого. Отметим, что информационная составляющая ФАК содержит всего 5 гармоник и объясняет более 97 % дисперсии ряда (табл. 2). На шумовые компоненты процесса приходится: на затухающие экспоненты - 0.7 % дисперсии, а на высокочастотные синусоиды -1.6 %.
Экстраполяционные возможности методики демонстрируются на рис. 5. Оценки параметров базисных функций были получены при значении лага Т = 30. Затем вычисленные параметры были использованы для прогноза ФАК на лаг, равный Т0= 60. Сравнение этих кривых показывает, что изображенное продолжение функции в основном отвечает интуитивным предположениям о ее возможной экстраполяции. На этом же рисунке изображена ФАК с Т0 = 60, полученная непосредственно по авторегрессионной модели. Сравнение полученных результатов показывает практически идентичное поведение функций на участке Т с [0,30]. Поведение ФАК, рассчитанной непосредственно по данным (ФАК1) и экстраполированной (ФАК2) на больших сдвигах существенно различно. Во-первых, ФАК1 по мере увеличения сдвига по сравнению с ФАК2 имеет тенденцию к уменьшению периода осцилляции, и во-вторых, наблюдается ее повышенное затухание. Все эти эффекты можно объяснить потерей точности расчета ФАК по фактическим данным при сдвигах, сравнимых с длиной обрабатываемой реализации. В то же время поведение экстраполированной ФАК определяется начальной, наиболее точной частью кривой и, следовательно, лишено этих недостатков.
к (!) !
4.2.2. Разложение ряда по базисным функциям
Для разложения последовательности по системе базисных функций воспользуемся алгоритмом (18)-(28). В качестве системы разложения примем некоторые из базисов, представленных в табл. 2, полученные при этом результаты изображены на рис. 6. Анализ показывает, что по мере увеличения числа компонент базиса аппроксимация довольно быстро приближается к исходному ряду, хорошо просматривается вклад отдельных гармоник базиса, а разложение по индивидуальным базисным функциям довольно хорошо аппроксимируются гармоническими колебаниями с переменными амплитудами.
Если принять разложение по базису 3-8 за информационную составляющую процесса, то можно получить шумовую компоненту наблюдений, изображенную на рис. 1. Из этого рисунка, в частности, видно, что помеха по форме имеет определенное сходство с информационной кривой с коэффициентом корреляции г 0.25. Шум, как и исходный процесс, имеет нестационарность по дисперсии.
Если воспользоваться рассмотренной выше аналогией между разностным уравнением и системой колебательных контуров, то полученную таким образом реализацию шума можно рассматривать, как возбуждающую линейную систему вынуждающую силу. Откликом на это внешнее возмущение является суперпозиция решений в виде затухающих экспонент или затухающих по экспоненте синусоид. Таким образом, анализ исходных данных по данной методике позволяет полностью восстановить модель процесса в непрерывном времени.
4.3. Гармонический анализ данных.
4.3.1. Гармонический анализ для случая априорно известных частот процесса
Возможности гармонического анализа данных используем для практически важного случая анализа сигналов, имеющих фазовую модуляцию. Обработка таких сигналов получила большое распространение в радиотехнике и системах измерений. При этом возникает целый ряд специфических проблем анализа и обработки данных, которые, на наш взгляд, наиболее эффективно могут быть решены в рамках разработанной выше теории. Известно, что фазово-модулированный сигнал в общем виде может быть представлен следующим выражением:
51(0 = Ассоъ[А ъхп(ат^-ф)\ (40)
где (р- мгновенная фазовая задержка; Ас - амплитуда гармонического сигнала; глубина модуляции; <ят -основная частота модуляции. Используя формулу косинуса суммы углов, выражение (40) может быть переписано в следующем виде:
(0 = АссоБ(^)соБ[^&\п(а>т^)]+АсБ\п(ф)Б\п[Р&\п(а>т^)\. (41)
Таким образом, фазовая задержка приводит к появлению сигналов, модулированных по косинусному и синусному законам. Обычно рассматривается следующая задача: по измерениям (40), представляющим сигнал ограниченной длительности Т, требуется определить фазу сигнала и оценить помехозащищенность измерений по отношению к шумам, действующим в измерительном тракте. Воспользуемся следующими известными соотношениями:
соб^Ш^)] = ^(р) + 2Ы2к(Р)со8(2кютГ)
к = 1,2,.,да, (42)
Бш[/Мп(®т0\ = 2Е/2к-1(^ш[(2к-1)®п/],
где Jn(P) - функции Бесселя первого рода порядка п и аргумента Д
Подставляя соотношения (42) в уравнение (41), получим следующую факторизацию сигнала (40), известную, как разложение Фурье:
51(0 = Ас[Зй(Р)со^(ф)+ 2 саь(ф)Ъ1:аШсо$(2кат£) +
2 яп(^2к-1(^т((2к-1)йт/)], к = 1,2,-,да. (43)
Выражение (43) показывает, что нечетные гармоники выходного сигнала зависят от а
четные - от соб(^>). При фиксированной глубине фазовой модуляции Р функции Бесселя превращаются в известные константы. Тогда выражение (43) будет представлять суммы гармонических колебаний с амплитудами, пропорциональными значениям косинуса или синуса фазы сигнала. Следовательно, для вычисления фазы достаточно определить амплитуды хотя бы одной пары четной и нечетной гармоник, присутствующих в сигнале. Очевидно, что точность определения фазы может быть улучшена, если использовать несколько пар таких гармоник с последующим осреднением результатов оценки. Таким образом, возникает типичная задача выявления гармонических колебаний по конечному набору данных.
Реализация описанного выше метода может быть достигнута как аппаратным путем, так и алгоритмически. В первом случае используется набор узкополосных фильтров, позволяющих выделять
гармоники сигнала, кратные основной частоте модуляции от,. Однако в процессе работы такого измерителя может наблюдаться аппаратный дрейф некоторых его параметров, ухудшающий точность измерений. Наиболее часто имеет место уход частоты и глубины модуляции от номинала, изменение в процессе работы параметров узкополосных фильтров, их неидентичность и т.д. Анализ показывает, что ошибки измерителя, связанные с дрейфом и неидентичностью его параметров, могут быть весьма существенными. Следовательно, возникает задача учета временного изменения параметров устройства с последующей их компенсацией. Учитывая, что фаза связана с непосредственно измеряемыми величинами существенно нелинейно, учет и компенсация дрейфа параметров становится сложной технической проблемой, намного удорожающей прибор. Другой важной проблемой при аппаратной реализации определения фазы является выбор времени, необходимого для осуществления единичного измерения Т. Этот параметр, с одной стороны, должен быть достаточно большим, чтобы смогли затухнуть переходные процессы узкополосных фильтров, а также уменьшить влияние разрывов ряда на краях интервала анализа. С другой стороны, время Т должно быть достаточно малым, чтобы уменьшить ошибки, связанные с естественным изменением измеряемого параметра за этот период.
Таблица 3. Значение функций Бесселя первого рода различных порядков для значений аргумента х = п и х = л/2
т 0 1 2 3 4 5 6 7
х = я -0.3042 0.2846 0.4854 0.3335 0.1514 0.05214 0.01455 0.003420
х = я/2 0.4720 0.5668 0.2497 0.06904 0.01400 0.002245 0.2983 10-3 0.3385-10-4
т 8 9 10 11 12 13 14
х = я 0.6961 10-3 0.1250 10-3 0.2009 10-4 0.2925-10-5 0.3891 10-6 0.4767-10-7 0.541310-8
х = я/2 0.3352-10-5 0.2946-10-6 0.2327-10-7 0.1669 10-8 0.1097 10-9 0.6649-10-11 0.3741 • 10-12
Указанных выше проблем можно избежать при алгоритмическом частотно-временном способе обработки измерений. Сущность предлагаемой методики заключается в следующем. Предполагается, что дрейф параметров измерителя представляет собой достаточно медленный процесс. Другими словами, полагаем, что этими изменениями можно пренебречь хотя бы на времени выполнения одного такта измерения. По гармоническому алгоритму определяются действующие в сигнале частоты. В идеальном случае, т.е. при отсутствии помех, эти частоты, согласно формуле (43), должны быть кратны частоте модуляции сигнала от. Таким образом, частота модуляции определяется на каждом цикле измерения и, следовательно, автоматически производится учет ее дрейфа. При наличии в измерительном тракте помех частота модуляции будет определяться с некоторой погрешностью. Ее уточнение может быть достигнуто за счет использования высших гармоник и учетом их кратности от,. Далее во временной области, в соответствии с (43), формируется система базисных функций, по которой выполняется разложение наблюденной функции. Коэффициенты полученного разложения используются для вычисления фазы сигнала. Фактически эта часть алгоритма идентична аппаратной узкополосной фильтрации сигнала. Однако при этом автоматически учитывается дрейф частоты модуляции, а фильтры имеют идеальные и идентичные характеристики. Следовательно, устраняется влияние тренда параметров фильтров и их неидентичность. Кроме того, как будет показано ниже, очень легко может быть учтен тренд параметра глубины модуляции. Таким образом, алгоритм определения фазы, основанный на гармоническом анализе, позволяет наиболее эффективно решать задачи обработки информации. Рассмотрим предлагаемую методику более подробно.
Для использования гармонического анализа необходимо представить сигнал (43) в дискретной форме. Пусть дискретизация выполнена через равные промежутки времени (А/). В этом случае будем иметь:
Sl(tk) = Лс[/о(Д)со8(ф)+ 2со$,(р)Ъ12к(Йсо$(2ка)т/к) + , ,„
к = 12 <х> (44) 2яп(^2ы(^т((2к-1)йт4)], """ '
где 4 = 0,1,...,Т - дискретный параметр времени; Т - время измерения функции, приходящееся на один цикл измерения фазы; = КА/.
Для определения периода дискретизации А/ воспользуемся теоремой Котельникова. Число измерений Т необходимо выбирать, исходя из периода самой низкой гармоники процесса (требуется, чтобы она имела в измерениях хотя бы один полный период).
Для практической реализации алгоритма определения фазы необходим расчет функций Бесселя первого рода. Эти функции могут быть вычислены через интегралы или бесконечные ряды. Наиболее просто на ЭВМ реализуется второй подход, который мы и рассмотрим более подробно. Функция Бесселя первого рода т-ого порядка с аргументом х выражается через ряд следующим образом:
Зт(х) = Ек[(-1)к(х/2)2к+т\/к!Г(т+к+1), к = 0,1,., да, (45)
где Г() - гамма-функция.
При реальных вычислениях бесконечная сумма (45) должна быть заменена на конечные пределы суммирования с оценкой достигнутой при этом погрешности. При х > 0 выражение (45) представляет знакочередующийся ряд и, следовательно, погрешность аппроксимации ряда конечной суммой может быть легко оценена (она равна первому отброшенному члену суммы). В табл. 3 приведены расчеты функций Бесселя нескольких порядков при значениях аргументов х = п и х = л/2. Для достижения точности £ = 0.00001 требуется удержать около десятка членов разложения (45).
Анализ этих таблиц, в частности, показывает, что функции Бесселя по модулю имеют резкий спад с увеличением их порядка. В связи с этим при глубине модуляции х = п может быть выделено около 7, а при х = л/2 - всего 5 гармонических компонент, имеющих более или менее значимую энергию. Все остальные гармоники имеют, по крайней мере, на два порядка меньшую энергию, чем основная энергонесущая частота процесса. Так как число выделенных гармоник в конечном счете определяет точность оценки фазы сигнала, глубину модуляции выгодней принимать равной х = п. Правда, с уменьшением параметра глубины модуляции растет (при одинаковых временах анализа) разрешающая способность метода. Отметим, что оценка числа значимых гармоник сигнала должна быть использована для определения периода дискретизации А/ по теореме Котельникова.
Выполним гармонический анализ данных и определим присутствующие в сигнале частоты. При этом порядок модели авторегрессии определим, исходя из числа энергетически значимых гармоник, присутствующих в процессе. Примеры фазово-модулированных сигналов, а также их комплексные и энергетические спектры в зависимости от величин фазовых задержек приведены на рис. 7-9. Анализ спектров показывает, что на частотах, кратных частоте модуляции, имеются хорошо выраженные спектральные пики. Максимумы в спектральной энергии, не связанные с этими частотами, объясняются конечностью выборки данных. Кроме того, из рисунков видно, что энергия процесса сосредоточена первых пяти гармониках, вклад высших частот в процесс сравним с энергией шума, обусловленного конечностью выборки. Уточним частотный состав процесса, используя кратности гармоник основной частоте от,. Анализ показывает, что этих данных достаточно для высокоточного определения гармоник модулирующей частоты в условиях малых помех. На основе полученных частот сформулируем следующую систему базисных функций:
ад,) = л^;
ад) = 2и7к(Р)со^(2ктт/д; к = 1,2,.,/, (46)
ад,) = 21/2Ы(^т[(2к-1)йт/1)\,
где / - число энергетически значимых гармоник ряда.
Легко видеть, что базис (46) имеет структуру, согласованную с разложением Фурье (44). При фиксированных а>т и Р этот базис имеет постоянную структуру, которая может быть запомнена, что оптимизирует время работы алгоритма. При дрейфе параметров базис должен периодически пересчитываться. Отметим, что может быть выбран и другой базис, отличный от (46). Например, в качестве базиса могут быть приняты любые две (или несколько) гармоник с частотами разной четности. Преимуществом базиса (46) по сравнению с другими является наиболее полный учет всей доступной информации и, следовательно, повышение точности вычислений. Однако самое важное преимущество этого базиса заключается в возможности учета тренда параметра /?. Разложив анализируемый ряд по базису (46), получим следующую факторизацию ряда (40)
ад) = Ас[А170(А)+А2ад,))+А3ад1)], (47)
где Аь А2, А3 - коэффициенты разложения функции по базису (46).
По последнему соотношению, учитывая выражение (43), легко определяем и оценку фазы сигнала
<р= шШ(А3/(А!+А2)). (48)
У(а>)
15-
Ю-
-35-1
1 ООО -
Рис. 7. Энергетические спектры фазово-модулированного сигнала при различных фазовых задержках: Ф=0 (1); ф=л/4 (2); ф=л/2(3) (Ща) - четная составляющая; У(а>) - нечетная составляющая спектра; £(®) - энергетический спектр)
Рис. 8.
а) ФМ-сигнал без помехи (ф=л/6); б) сигнал с адитивной помехой (Аш=0.5); в) результат частотной фильтрации
Рис. 9.
а) ФМ-сигнал без помехи (ф=л/6); б) сигнал с адитивной помехой (4ш=1); в) результат частотной фильтрации
Рассмотрим результаты моделирования по указанной выше теории. При моделировании фаза выбиралась случайным образом из интервала [0; 2л], а значения сигнала (при фиксированных от и /3) вычислялись по формуле (40). Некоторые реализации функции при различных значениях фазы (/3 = п, Ас = 1, от = 0.25) представлены на рис. 8, 9. Полученные данные подвергались затем гармоническому анализу для определения частоты и фазы. Проанализируем полученные результаты.
Точность оценки фазы не зависит явно от числа учитываемых при анализе гармоник, т.е. не определяется полнотой базиса (46). В то же время экспериментами установлено, что в идеальном случае (при отсутствии помех) описанный выше алгоритм позволяет определить фазу сигнала в пределах точности арифметики ЭВМ. При этом длина выборки данных должна быть не менее чем два периода основной гармоники при восьми отсчетах на период на самой высокой частоте учитываемой синусоиды. Базис (46) включал в себя информацию по семи первым гармоникам. Точность аппроксимации функции также является высокой (коэффициент корреляции между сигналом и его приближением (47) составлял более 0.995). При использовании неполных базисов высокая точность оценки фазы также сохраняется, однако ухудшается точность аппроксимации.
Этот результат является вполне объяснимым для рассматриваемого идеализированного случая: все вычисления детерминированы и однозначны и допускают ошибку только за счет машинной точности ЭВМ. Важным является и то, что легко поддаются оценке ошибки, связанные с представлением бесконечного ряда (45) конечной суммой. При этом точность оценки фазы в идеальном случае не зависит в конечном счете от точности аппроксимации ряда.
Рассмотрим возможности гармонического анализа при присутствии в данных помех и вариации параметра глубины модуляции. Будем по-прежнему считать, что число гармоник процесса известно (т.е.
известен оптимальный порядок регрессионной модели). В этом случае модель данных с учетом возможных шумов в измерительном канале будет иметь следующий вид:
£(4) = А^ъ№т(аЛк-ф)\+Ашг1Цк),
(49)
где п(/,) - шум наблюдений в дискретных точках единичной амплитуды; Аш - амплитуда шума; /к -параметр дискретного времени.
Модельный шум единичной амплитуды генерировался согласно следующему выражению:
п(4) = иг(4)-0.5,
(50)
где иг(4) - случайное число с равномерной плотностью распределения в интервале [0; 1\.
Для реализации временной последовательности значений шума использовался стандартный генератор псевдослучайных чисел, входящий в программное обеспечение ЭВМ. Типичная реализация используемого при моделировании шума приведена на рис. 10, а его функция автокорреляции и спектр мощности изображены, соответственно, на рис. 11, 12. Анализ ФАК шума показывает, что она имеет значимый отсчет только на единичном сдвиге, а на более высоких значениях лага имеет относительно небольшие величины (зависимость типа дельта-функции). Такое поведение ФАК свидетельствует о том, что используемый генератор случайных чисел действительно генерирует сигнал, близкий по свойствам к белому шуму. Об этом же говорит и спектр мощности сигнала (рис. 12). Вектор шума складывался затем с отсчетами сигнала. Считается, что полученный таким образом временной ряд и представляет модель наблюдений (49).
Рассмотрим влияние шума наблюдений на работу алгоритма определения фазы сигнала. Для этого генерировался шум наблюдений по соотношению (50) при различных значениях параметра Аш. Добавляя полученный шум к идеальному сигналу (40), имеем зашумленный сигнал (49). Далее применим к полученному сигналу алгоритм определения фазы и оценим полученную погрешность оценки фазы сигнала.
Рис. 10. "Типичная" реализация шума с равномерной функцией распределения (4ш=1)
Рис. 11. ФАК реализации шума
Рис. 12. Спектр мощности реализации шума
На рис. 9 приведен сигнал с помехой, равной Аш = 1 (^= л/3). Ошибка вычисления фазы составила в этом случае 8(р = 0.0136 (рад.). Рис. 8 демонстрирует сигнал и сигнал с помехой с уровнем шумов Аш = 0.5. Задавая разные уровни модельного шума и фазовые углы, получим соответствующую
статистику. Моделированием установлено, что абсолютная ошибка определения фазы за счет действия белого шума не зависит прямым образом от значения измеряемого параметра. В то же время с ростом помех пропорционально растет и ошибка определения фазы. Полученная статистика затем была обработана по методу наименьших квадратов (МНК) для получения коэффициентов линейной регрессии. Для принятых условий связь ошибки определения параметра с величиной аддитивной помехи имеет следующий вид:
8(р= 0.0136 Аш. (51)
Эта зависимость была промоделирована для амплитуды шума в интервале Аш с. 0^2, для различных значений измеряемой фазы и реализаций помехи. Установлено, что связь ошибки в оценке фазы с действующими в измерительном канале помехами имеет устойчивый и линейный характер. Таким образом, наличие в измерительном канале белого шума приводит к систематическому смещению (завышению) оценки фазы сигнала. Как было показано выше, оценка параметров шума легко может быть получена в рамках предлагаемой методики и, следовательно, оценка фазы в этом случае легко может быть уточнена по соотношению (51).
Следующая группа экспериментов заключалась в определении влияния вариации глубины модуляции Р на точность определения фазы сигнала. Обозначим номинальное значение глубины модуляции через Д>, а ее отклонение от номинала - 8f3. Таким образом, принимаем, что текущее значение глубины модуляции Нравно f3= f30+8f3. Выполним разложение функции Бесселя в ряд Тейлора и, ограничившись первыми двумя слагаемыми разложения, будем иметь
JÁP) = Jn(fr) + ЛШ 8р. (52)
Оценка производных функций Бесселя может быть выполнена численно или аналитически. В последнем случае для этого можно воспользоваться известными рекурентными соотношениями
Jn (x) = Jn-1 (x) - (n/x)Jn(x). (53)
Подставляя выражение (52) в (43), получим следующее разложение сигнала по гармоникам: Si(ti) = Ac[J0($,)cos(^) + 2cos( <p)'EJ2k(A))cos(2k&mti) + 2sin(^) J_i($,)sin((2k-1)®„/,)]+
+ Ac[J'0(P0)cos(^)+2cos(^)U2k'(P0)cos(2k&mt¡)+ 2sin(p)I/2W^sin((2M)fimt)]¿A (54)
Первое слагаемое выражения (54) является значением сигнала при номинальной величине глубины модуляции. Если разлагать сигнал (49) по базису (46), то будет выделена именно эта компонента данных. Второе слагаемое представляет вклад в ряд, обусловленный уходом параметра от номинала. Следовательно, дрейф параметра fí приводит к неточному восстановлению функции по базису (46), причем ошибка рассогласования пропорциональна отклонению 8f3. Для того, чтобы учесть тренд параметра ¡3, сформируем новую систему базисных функций следующего вида
ад) = J0(A); Y1(t) = 2 U2k(^)cos(2kümt,); ym = 2 J ($,) sin[(2k-1)©„/,];
Ys(ti) = 2 L/2k'(^)cos(2k©mt,); Y4(t,) = 2 Jw'^) sin[(2k-1)©„/,].
Структура этого базиса становится ясной после его сравнения с разложением (54). Разложив сигнал по данному базису, получим следующее его представление:
S1 (t) = AcA Y()(ti)+A2Y\ (tl)+A3Y2(tl)+A4Y3(tl)+AiY4(tl)], (56)
где Ai (i = 1,2,...,5) - коэффициенты разложения сигнала по базису (54).
В идеале коэффициенты этого разложения будут соответственно равны: 1) cos(^>); 2) cos(^>); 3) sin(^>); 4) cos(p)Sj3; 5) sin(^»)¿¡/?. Эти коэффициенты могут служить для раздельной оценки фазы, а также для расчета отклонения глубины модуляции от номинала. Анализ соотношений (55), (56) показывает, что оценка фазы сигнала не зависит прямым образом от изменения параметра глубины модуляции и может рассчитываться по первым трем коэффициентам разложения. После оценки глубины модуляции по четвертому и пятому коэффициентам разложения может быть рассчитано ее отклонение от номинала. Если ограничиться в разложении (56) только первыми тремя членами, то точность приближения функции будет прямым образом зависеть от параметра 8f3. Следовательно, эта невязка является характеристикой вариации глубины модуляции (а также действия в измерительном тракте шумов).
Хотя рассмотренная выше модель не устанавливает зависимости между ошибкой в оценке фазы сигнала и дрейфом параметра 8f3, экспериментами установлено, что такая связь существует. Эта
зависимость обусловлена конечностью анализируемой выборки данных, в результате чего начинают проявляться так называемые граничные эффекты. Поэтому качество разложения по базису (55) зависит от поведения функции на концах интервала анализа. Ошибки, связанные с краями, будут малы, если на концах интервала не будет наблюдаться слишком больших разрывов функции. Отметим, что аналогичные эффекты имеют место и при Фурье-анализе (так называемое явление Гиббса). Следовательно, для уменьшения такого рода ошибок желательно, чтобы период наблюдения функции был бы кратен частоте модуляции. Однако строго выдержать это требование по реальным данным затруднительно.
Рассмотрим задачу улучшения оценки фазы при изменении глубины модуляции несколько с другой стороны. Чувствительность ошибки в оценке фазы, в зависимости от вариации глубины модуляции, связана со свойствами функций Бесселя в районе номинального значения параметра и длины анализируемой реализации. Аналитическое исследование этих зависимостей, с учетом их существенной нелинейности, является трудной задачей. Поэтому проведем этот анализ численными методами.
Будем варьировать параметр глубины модуляции сигнала с одновременным определением ошибки в оценке фазы сигнала, вызванной этим трендом. Накопленная таким образом статистика была обработана по МНК с целью определения регрессионной связи этих величин. В результате установлено, что для Д) = п и выбранного периода наблюдений (рис. 8, 9), ошибка оценки фазы имеет следующую связь с изменениями глубины модуляции
д(р = 0.033 др. (57)
Как видно из выражения (57), ошибка в определении фазы из-за ухода глубины модуляции от номинала является линейной функцией. Следовательно, измерительная схема является устойчивой по отношению к дрейфу этого параметра. Улучшить оценку фазы можно непосредственно по уравнению связи (57). Следует, однако, отметить, что при других соотношениях между периодом наблюдений и периодом основной гармоники модуляции, численное значение коэффициента регрессии (57) будет иным. Поэтому этот подход не является достаточно универсальным. Отметим также, что эта же проблема существует и при аппаратной реализации алгоритма.
Некоторое улучшение оценки фазы может быть достигнуто за счет удержания в разложении Тейлора (52) квадратичной (и выше) компонент. Экспериментами, однако, установлено, что вклад квадратичных членов слишком мал и не улучшает существенным образом оценку фазы. Другой способ улучшения оценки заключается в использовании в алгоритме не априорно заданного параметра глубины модуляции, а его фактического значения, полученного в результате той или иной его оценки. Этот подход может быть выполнен в итерационном цикле. На первом этапе за основу принимается номинальное значение этого параметра, полученное по априорной информации или по результатам предыдущих измерений. Затем это значение уточняется за счет вычисления невязки др по алгоритму (5256), по этому уточненному значению производится новый пересчет данных и т.д. Критерием достижения нужного результата может служить, например, невязка аппроксимации функции по (56).
Следующим способом улучшения оценки фазы является учет ошибки за счет использования зависимости (57). Для этого необходимо выполнить моделирование для установления зависимости коэффициента регрессии от соотношения периода наблюдений и глубины модуляции. Во всяком случае, увеличение точности оценки фазы сигнала при изменении параметра Р не представляет принципиальных трудностей. Отметим также, что выражение (54) можно использовать в качестве алгоритма для определения знака и модуля отклонения параметра модуляции от номинала. Этот алгоритм может быть использован для аппаратной коррекции параметра модуляции.
5. Заключение
Таким образом, применение для анализа исходных данных теории линейных динамических систем дает возможность эффективного решения многих актуальных задач обработки информации. Получаемые при этом результаты имеют ясный физический смысл и легко поддаются интерпретации. Разработанная выше методология обработки одномерных временных данных легко может быть обобщена на векторные временные ряды и многомерные поля.
Литература
Высоцкий В.М., Ивахненко О.Г., Чеберкус В.И. Долгосрочное прогнозирование колебательных процессов с помощью выделения гармонического тренда оптимальной сложности по критерию баланса переменных. Автоматика, № 1, с.29-34, 1975.
Губанов B.C. Обобщенный метод наименьших квадратов (теория и применение в астронометрии). СПб., Наука, 319 е., 1997.
Драница Ю.П. Об одном методе моделирования нестационарных динамических систем и процессов. Вестник МГТУ, т. 3, № 1, с.55-66, 2000.
Корн Г., Корн Т. Справочник по математике (для научных работников и инженеров). М., Наука, 832 е., 1977.
Котюк А.Ю., Ольшевский В.В., Цветков Э.И. Методы и аппаратура для анализа характеристик случайных процессов. М., Энергия, 237 е., 1967.
Pao С.Р. Линейные статистические методы и их применение. М., Наука, 73 е., 1968.
Ткачев Г.Н., Ланкис Л.К. Экспериментальное опробование методики пересчета вариаций естественного электромагнитного поля Земли по функциям электропроводности. Геофиз. журнал, т. 15, № 4, c.76-80, 1993.
Троян В.М., Соколов Ю.М. Методы аппроксимации геофизических данных на ЭВМ. Л., Ленинградский университет, 302 е., 1989.
