Спектральный анализ на основе линейной динамической модели

УДК 517.958+517.962.2

Спектральный анализ на основе линейной динамической модели

12 2 Ю.П. Драница , А.Ю. Драница , О.В. Алексеевская

1 Политехнический факультет МГТУ, кафедра высшей математики

и программного обеспечения электронно-вычислительных машин

2

ЗАО "Ланит ", г. Москва

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

Abstract. The task of parametric estimation of the spectral density with a high resolution on the basis of the linear dynamic model has been put and solved in the paper. The high resolution of the method has been obtained due to extrapolation of the process autocorrelation function throughout the number axis. The estimations of necessary dynamic parameters of the process (frequency and factors of fading) have been performed using the current measured data. The received spectral estimations are being expressed by simple analytical formulae and allow to interpret the results of spectral estimations at meaningful level. The methodics gives an opportunity to obtain individual spectral estimations automatically for any frequency, this is rather useful for more delicate research of the process structure.

Ключевые слова: линейная модель, спектр, высокое разрешение, аналитическое представление спектра, экстраполяция, функция автокорреляции

Key words: linear model, spectrum, high resolution, analytic representation of spectrum, autocorrelation function

1. Введение

Одним из способов изучения временных рядов является вычисление спектральной плотности, которая показывает распределение энергии процесса по частотам колебаний. Отметим, что в классическом математическом анализе спектральный подход практически не используется, а в математической статистике предпочтение отдается временным рядам.

Роль спектральных методов возросла после того, как были выяснены причины статистической несостоятельности спектральных оценок (Теребиж, 2005) и предложены процедуры, решающие эту проблему. Однако роль спектральных оценок при научных исследованиях по-прежнему незаслуженно мала. По мнению авторов, важность спектральных методов заключается в том, что спектральный состав процесса тесным образом связан с его функцией автокорреляции (ФАК). Эта связь формально выражается теоремой Винера и Хинчина (Краус, Вошни, 1975) и позволяет по ФАК вычислить спектральную мощность и наоборот. В свою очередь, ФАК является основополагающим понятием теории линейных систем, позволяющей синтезировать динамическую систему. В связи с изложенным, задача адекватной оценки функции спектральной плотности является актуальной как с научной, так и практической точек зрения. Решению данной задачи с позиций динамического линейного моделирования и посвящена данная работа.

2. Общая постановка проблемы

Пусть имеются измерения временной последовательности некоторых физических величин, выполненные в n точках, с постоянной дискретностью At. Эти измерения образуют временной ряд

X =

Для упрощения дальнейших выкладок будем считать, что ряд (1) центрирован. Отличие временных рядов от обычной совокупности n измерений заключается в упорядоченности отсчетов временного ряда. Это значит, что для временного ряда нужно задавать две связанные между собой последовательности: отсчеты

времени и соответствующие им значения измеряемой величины. Обычно в качестве упорядочивающего фактора выступает время, но могут быть другие одномерные или многомерные переменные.

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

Будем исходить с точки зрения статистического подхода и трактовать данные (1) как временную выборку некоторой случайной величины. Если рассматриваемый временной ряд обладает инвариантностью относительно основных статистик второго порядка, то на выборке (1) можно определить эмпирическую оценку r функции автокорреляции (ФАК) по следующей формуле

rТ) = -ТI , т = ±1, ±2,.. ±1,

где l - максимальная задержка (или лаг) сигнала, используемая при оценке ФАК. Ввиду свойств симметрии ФАК (2) является четной функцией: г(т) = г(-т) = rT, поэтому для ее вычисления достаточно ограничиться нахождением г(т) при т > 0. В статистике показано, что оценки (2) являются несмещенными и состоятельными.

Одним из методов линейной обработки стационарных временных рядов является спектральный анализ, позволяющий охарактеризовать частотный состав измеряемого сигнала, который показывает распределение энергии по частотам колебаний - так называемую спектральную плотность, или спектр мощности. Преобразование Фурье является математической основой, связывающей временной сигнал с его представлением в частотной области. Как показано в (Теребиж, 2005), в связи с ограниченностью выборки (1) и присутствием в данных шума измерений спектральные оценки относятся к так называемым обратным и некорректно поставленным задачам.

Результаты спектральных оценок существенным образом зависят от принятой модели данных. В соответствии с этим принято различать непараметрические и параметрические спектральные оценки. В первом случае для оценок используется только последовательность измерений (1). Во втором случае используются некоторые априорные предположения, например, предполагается, что для данных (1) может быть использована авторегрессионная модель p-ого порядка (краткое обозначение - AR(p)). Отметим, что оценка спектральной плотности на основе AR(p) модели сводится к оценке p+1 неизвестных параметров (Теребиж, 2005).

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

AtAm > 2п.

Допустим, что время регистрации сигнала составляет T. Тогда, в соответствии с принципом неопределенности (3), разрешающая способность (Bp) классических оценок спектральной плотности оценивается величиной Bp ~ 2n/T, т.е. обратно пропорциональна времени наблюдения. На непараметрическом оценивании базируется классический спектральный анализ, основанный на преобразовании Фурье. Для улучшения качества спектральных оценок в классическом спектральном анализе прибегают или к их усреднению или сглаживанию с помощью различных окон (Хемминг, 1987).

Более высокое разрешение дают параметрические методы оценки спектральной плотности (Марпл, 1990; Теребиж, 2005). Стабилизация этих оценок обеспечивается за счет привлечения дополнительной априорной информации, например, предполагается, что данные (1) могут аппроксимироваться авторегрессионной моделью. Спектральная плотностьf(ffl|a) на основе AR(p) модели имеет вид

_ 2

f

11 -1am exp(-'2ка>ш) |

где а^ - дисперсия шума, порождающего авторегрессионный процесс, вектор а - коэффициенты АЯ(р) модели, а временной шаг принят равным 1. Ясно, что оценить всего р+1 неизвестных параметров, полностью задающих функцию (4), гораздо легче, чем все значения спектральной плотности на достаточно плотной сетке частот. Рассмотрим, почему параметрические методы имеют более высокое разрешение, чем непараметрические.

В классическом подходе считается, что за пределами интервала наблюдения исследуемая зависимость равна нулю. Отметим, что дополнение данных нулями вводит искусственный разрыв в аппроксимируемую функцию, который приводит к сдвигу спектра в сторону высоких частот. Если имеющийся сегмент данных короче ФАК анализируемого стационарного процесса, параметрические методы производят ее экстраполяцию за предел эмпирической оценки ФАК, что и приводит к увеличению разрешающей способности спектральных оценок. Таким образом, разрешающая способность спектральных оценок параметрическими методами зависит от возможности адекватной экстраполяции ФАК.

Существующие методы параметрического спектрального оценивания, основанных на AR(p) моделях, базируются на формуле (4). Нами предлагается спектральное оценивание, основанное на других принципах, а именно на следующем. На основе линейной модели производится оценка динамических параметров процесса (частоты, коэффициенты затухания). Эта информация служит основой для создания базиса разложения, в котором сначала производится аналитическая аппроксимация ФАК, а затем ее экстраполяция на временной интервал (-<»,<»). Полученное продолжение ФАК далее используется для аналитической оценки спектральной плотности. Получаемые в результате спектральные оценки выражаются суммой простых аналитических выражений. Это, по нашему мнению, придает спектральным оценкам большее физическое обоснование, чем вычисления по формуле (4), и позволяет выполнить качественную интерпретацию спектральных оценок. Кроме того, предлагаемый метод позволяет автоматически получать индивидуальные спектральные оценки для любой частоты, что является весьма полезным при более тонком изучении процесса.

3. Предлагаемое решение проблемы

3.1. Задача аналитической аппроксимации ФАК

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

Y1 (0 = ехра 0,..., Yj = ехр(аД

Yj+1

Yj+2m-l

где j - число базисных функций в виде затухающих экспонент; l - общее число функций в базисе; t -дискретный параметр времени, ш - общее число ненулевых частот. Отметим, что каждой ненулевой частоте соответствует две функции вида Yk(í) = ехр(ак0 соб^О, Yk+1(t) = ехр(ак0 бш^О. Для удобства применения все базисные функции упорядочены по возрастанию частоты.

Сформированный таким образом базис позволяет достаточно эффективно решать ряд прикладных задач таких, например, как аппроксимация и экстраполяция функций, оценка шумовой компоненты сигнала, его разложение по системе функций и другие задачи. Мы будем рассматривать задачу разложения ФАК в этом базисе. Введем вектор Я=(г0, гь..., гТ)Т, где (*)г - операция транспонирования, и матрицу значений базисных функций следующего вида

А =

^

^

^

Yl

где ш - общее число функций, которые будут использоваться при разложении; Т - максимальная задержка (лаг) ФАК. Разложение ФАК по этому базису в матричных обозначениях будет иметь следующий вид

Я = А^

где Z - вектор столбец коэффициентов разложения (весовые коэффициенты). Сумму квадратов невязки между аппроксимацией (6) и эмпирической оценки ФАК можно представить следующим образом

=

Проблема состоит в определении оптимальных (по некоторому критерию) значений вектора Z. При ш < T аппроксимацию (6) можно рассматривать как типичную задачу метода наименьших квадратов (МНК) с критерием оптимальности в виде минимума суммы квадрата невязки (7). Нам необходимо найти такие оценки вектора Z параметров модели (6), которые удовлетворяли бы этому условию. Для

раскрытия экстремума функций Б(г) воспользуемся правилами векторного дифференцирования (Рао, 1968) и найдем необходимое условие экстремума

дБ/дг = - К Я + Ялг = 0,

откуда имеем

Алг = К А.

Условием минимума квадратичной формы £ означает, что матрица

д-БШ2 = А л

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

г =

где (• )-1 - операция вычисления обратной матрицы. Полученное решение (9) позволяет представить ФАК следующим образом

г

где коэффициенты А,, В, - это составляющие вектора г. Соотношение (10) можно представить в несколько ином виде

г

где ф, = аг^(В,/А,) - начальная фаза компоненты разложения; Сi = (А,2 + В,2)05 - обобщенная амплитуда. Отметим, что параметр 70 может удовлетворять как условию 70 < 7, так и 70 > 7. В первом случае производится интерполяция сигнала, а во втором - его экстраполяция или прогноз. Из выражения (10) следует, что компоненты вектора Z представляют собой весовые коэффициенты, с которыми базисные функции входят в разложение ФАК. Для нулевых частот (щ = 0) коэффициенты В, в формуле (10) равны нулю и, следовательно, начальные фазы в разложении (10) также могут быть приняты (ф, = 0). Следовательно, нулевые частоты являются частным случаем косинусного представления и входят в разложение (11). Очевидно, что данная аппроксимация может быть использована и для любой другой функции ф(к), имеющей ровно 70 отсчетов. Для этого достаточно подставить эту функцию в уравнение (6) вместо вектора Я и получить решения (10) или (11).

При анализе коротких рядов данных часто возникает задача надежной оценки корреляционных функций. Известно, что достоверная оценка ФАК возможна только для ее смещений в 5-10 раз меньших, чем длина анализируемого ряда (Котюк и др., 1967). При больших сдвигах оценки становятся неустойчивыми и могут не удовлетворять нужным свойствам. Одним из важнейших требований, которому должна удовлетворять эмпирическая оценка ФАК, является ее положительная определенность (Теребиж, 2005).

Симметричная матрица л называется положительно определенной, если для всех векторов х справедливо неравенство х7лх > 0. Только в этом случае ФАК позволяет адекватно рассчитывать спектр мощности и может быть использована в процедуре МНК. Нарушение положительной определенности матрицы сигнализирует о ее плохой обусловленности и численной неустойчивости обращения. В этом случае, например, амплитудный спектр может иметь отрицательные компоненты, а оценки параметров линейной модели не будут устойчивыми и обладать оптимальными свойствами. Другими словами, такой ФАК может не соответствовать реальный выход линейной системы. Отметим, что аппроксимация ФАК по базису (5) автоматически гарантирует положительную определенность этой функции.

3.2. Расчет спектральной плотности компонент разложения

Как известно, спектральная плотность Б(щ) и корреляционная функция г(т) связаны между собой парой преобразований Фурье:

г

Б

где , - мнимая единица. Формулы (12, 13) представляют собой комплексную форму прямого и обратного преобразования Фурье. Расчет спектральной плотности по формуле (13) для эмпирически определенной

функции корреляции обычно выполняют численным интегрированием. При факторизации ФАК по базису затухающих экспонент и гармоник появляется возможность вычисления спектра (13) аналитически. Допустим, что выполнена аппроксимация ФАК суммой затухающих экспонент и затухающих по экспоненте гармоник (10). Рассчитаем спектр составляющих типа затухающих экспонент, ФАК которых представляется выражением

г(т) = Аке

-щ М

где Ак, ак - соответственно амплитуда и коэффициент затухания экспоненты. Подставляя выражение (14) в формулу (13), будем иметь:

Бек(ю) = А |е-ак= А |а+1е-акТв-мс1т ]=

" 0 ^ |еак —1ю)Т ¿Т +1е<ак +1ю)Т¿Т ] .

=А

2п

Интегралы (15) вычисляется элементарно

5ек (ю) = А 2%

А

2п

1

(ак —1ю)Т -(ак +1ю)Т

е Iо е | --- — (-

ак—1ю " ак+1ю

= Акак /(к(ак + ю2).

- +-

ак — 1ю ак + 1ю

к

Спектр (16) представляет четную одномодальную функцию с максимумом на нулевой частоте. Спад функции к нулю при \юк\ ^<» определяется коэффициентом затухания ак. При уменьшении ак в спектре большой вес начинают приобретать низкие частоты. При увеличении ак веса высоких частот начинают увеличиваться, при ак ^ <» спектр приближается к равномерному (так называемому "белому") спектру, в котором нет преобладания каких-либо частот.

Рассмотрим спектр процесса, имеющего ФАК в виде затухающего по экспоненте косинуса или

синуса:

г

г(т) = Аке-щщ бШ^Т). Косинус и синус в комплексной форме могут быть представлены с помощью формул Эйлера:

соБ(юкт) = (ею + е-ШкТ )/2, Б1п(юк т) = (ею — е)/21.

Для расчета спектра затухающего косинуса, подставим первую формулу выражения (17) в формулу (13). В результате, разбив пределы интегрирования на отрицательное и положительное время, будем иметь два несобственных интеграла следующего вида:

5»(ю)=А

| еакТ (е'юкТ + е~'щТ )е~ш ¿Т +1 еа (е1ЩТ + е~,щТ )е~ш ¿Т

Рассмотрим интеграл для отрицательного времени

—,юкТ

1с

0 0 „аТ I ЛюъТ , _ —1юк

(ю) = | еакТ (е'юкТ + е "ЩТ )е~гютёт = | еатеюе~гютёт + | еате'10^е~1югёт =

00 = | е(ак +1(юк—ю))Т ¿Т + | е(ак —1(юк +ю))Т ¿Т

который вычисляется элементарно

1сп(ю) = ехр(ак + 1(юк —ю) |0 + ехр(ак — 1(юк + ю) |0 =

ак + 1

ак — 1

1

1

= ак — 1

ак + 1(юк — ю) ак — ¡(юк +ю) ак 2 + (юк — ю)2 ак 2 + (юк +ю)2

0

Аналогичным образом рассмотрим интеграл для положительных времен

ию) = |е~акТ (еЩкТ + е~ШкТ)е-щ1ёт = ]еакте,щте-тёт + ]е~акТе~ШкТе~гатёт -

= | е( ~ак + Щк ~Щ))ТёТ+| е( ~ак - (Щк +Щ))ТёТ

который вычисляется элементарно

1сп(ю) = ехр(~ак + ,(щк -щ) I» + ехр(~ак - ,(щк + щ) I» =

-ак + 1(ак -щ) -ак -+ щ) '0

= -1 -1 = ак + ¡

-ак + 1(юк -Щ) -ак - 1(юк + Щ) ак2 + (щ -щ)2 ак2 + (щ + щ)2

Подставив в выражение (19) значения интегралов (20) и (21), после элементарных преобразований получим следующее выражение спектральной плотности затухающего по экспоненте косинуса:

Бек Щ) = А-

2п ак + (щ -щ)2 ак + (щ + юу Зависимость (22) представляет двухмодальную четную функцию с экстремумами на частотах ±щк. Вид графика (22) зависит от соотношения параметров ак, щ, т. е. оттого, что преобладает в ФАК: убывание по экспоненте или гармоническое колебание. При малом ак процесс практически ведет себя как гармоническое колебание, и спектр представляет собой сравнительно узкую линию на частоте этой гармоники. При увеличении ак спектральная линия начинает размываться, становясь все более широкой. В пределе при ак^» спектр сигнала приближается к "белому" спектру.

Для вычисления спектра процесса, имеющего ФАК вида затухающего по экспоненте синуса, воспользуемся второй формулой выражения (18). Подставив эту формулу в (13), будем иметь

(щ_А

0

| еакт (ещкт - е-щт )е~'стёт + ] е а (ест - е~щт )е-шёт

-» 0

Аналогично предыдущему случаю, рассмотрим интеграл для отрицательного времени

0 0 0

иш) = | еакТ (егс'кТ - е~гс'кТ )е-ютёт = | еакТеЩкТе-ютёт - | еакТе~ШкТе-ютёт -

0 0 = | е(ак +,(к-Щ))тёт- | е(ак-,(к +Щ))тёт

который вычисляется элементарно

1сп

= 1 1 = ак - ,

ак + 1

-ак + - щ -ак - + щ ак + (щ - щ)2 ак + (щ + щ)2 Проведем следующие преобразования

_ Т _~Шк +

П 2 , \ 2 2 , \ 2

, ак +

_ _-,ак - (щ -щ) -,ак + (щ +ю).

^р 2 2 2 2

,' ак +

В результате выражение (23), с учетом соотношений (26), будет иметь следующее выражение спектральной плотности затухающего по экспоненте синуса:

^ О) = А-

2п ак +

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

5к

которое показывает распределение энергии затухающей гармоники по частотам. Спектр затухающей гармоники простирается от -да до да с резонансной частотой ±юк. Выражения (16, 22) представляют собой четную, а (27) - нечетную составляющие комплексного спектра. Суммарные четные и нечетные составляющие комплексного спектра ФАК, вследствие принципа линейности, будут, соответственно, равны

Б

кйл

щ) = ^ {-г^1-г +-),

к~12л ак + (щ -щ)2 ак + (щ + щ)2

где тп - число ненулевых частот. Отметим, что в выражении (30) и во второй сумме (29), индекс суммирования к перебирает последовательность ненулевых частот. Выражения (29) и (30) позволяют вычислить спектральную плотность процесса в следующем виде

Бк(ю) = Бч(т)2 + Бн(ю)2,

которое показывает суммарное распределение энергии сигнала по различным частотам. Отметим, что полученные оценки спектральной плотности имеют четкий физический смысл, обусловленный представлением данных (1) выходом некоторой линейной системы. Повышение разрешения спектральных оценок достигнуто за счет применения линейной модели и экстраполяции компонент ФАК на бесконечный временной интервал.

4. Заключение

Спектр процесса можно непосредственно вычислить, применив формулу (13) к исходной временной последовательности. В результате будет получена так называемая периодограмма Шустера. Бартлет (Теребиж, 2005) показал, что для периодограммы Шустера расхождение между истинной спектральной плотностью и ее оценкой не уменьшается с ростом объема выборки п, т.е. оценки являются несостоятельными, а сама постановка проблемы - некорректно поставленной задачей. Поискам состоятельных и эффективных оценок спектра мощности были посвящены усилия многих специалистов в течение столетия после появления пионерской работы Шустера.

Одним из способов построения состоятельных оценок спектра мощности является их сглаживание - классический подход в случае некорректных постановок. С другой стороны, в статистике показано, что эмпирические оценки ФАК являются несмещенными, эффективными и состоятельными. Более того, само вычисление ФАК является сглаживающей операцией. Следовательно, все предпосылки для получения состоятельных оценок выполнены, и спектр мощности, вычисленный по ФАК, должен обладать этими же свойствами.

Другое дело, что дисперсия оценок ФАК по мере увеличения задержек сигнала возрастает вследствие уменьшения объема выборки. Решение этой проблемы, в рамках предлагаемого подхода, заключается в следующем. Оценка ФАК по экспериментальным данным выполняется для таких задержек сигнала, при которых ее расчеты являются еще статистически обоснованными. По этим относительно точным оценкам выполняется аппроксимация ФАК и ее экстраполяция с последующим вычислением спектра мощности. Такой подход, по мнению авторов, позволяет на объективном уровне решить все проблемы адекватного спектрального анализа.

Литература

Драница Ю.П., Драница А.Ю. Корректный метод оценки спектральной плотности на основе линейной модели. Матер. Междунар. научно-технич. конф. "Наука и образование-2007". НТЦ "Информрегистр" 0320700491 от 05.03.2007, Мурманск, МГТУ, с.116-120, 2007.

Драница Ю.П., Драница А.Ю. Некоторые аспекты интерпретации экспериментальных данных на

основе теории линейных динамических систем. ВестникМГТУ, т.12, № 1, с.125-131, 2009. Клаербоут Ф. Теоретические основы обработки геофизической информации. М., Недра, 304 с., 1981. Краус М., Вошни Э. Измерительные информационные системы. М., Мир, 310 с., 1975. Котюк А.Ю., Ольшевский В.В., Цветков Э.И. Методы и аппаратура для анализа характеристик

случайных процессов. М., Энергия, 237 с., 1967. Марпл-мл. С.Л. Цифровой спектральный анализ и его приложения. М., Мир, 584 с., 1990. Рао С.Р. Линейные статистические методы и их применение. М., Наука, 73 с., 1968. Теребиж В.Ю. Введение в статистическую теорию обратных задач. М., Физматлит, 376 с., 2005. Хемминг Р.В. Цифровые фильтры. М., Недра, 222 с., 1987.