П2.10 Параметрический спектральный анализ

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

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

Наибольшее распространение получил подход, основанный на построении линейных стационарных моделей. В соответствии с этим подходом анализируемый процесс рассматривается как процесс на выходе линейной системы, на вход которой поступает белый шум. Поскольку белый шум имеет равномерную спектральную плотность, то, пропуская его через линейную систему бесконечной размерности (или, что то же самое, через бесконечное число линейных систем, соединенных последовательно или параллельно), можно получить спектральную плотность мощности любой формы. Поэтому линейная модель пригодна для аппроксимации многих встречающихся на практике детерминированных и случайных процессов с дискретным временем. Спектральная характеристика линейных моделей полностью определяется параметрами модели и дисперсией белого шума $\sigma_{\xi}^2$.

Весь класс линейных моделей может быть представлен в виде:

$$x(i)= - \sum_{k=1}^{p}a_k x(i-k) + \sum_{k=0}^{q} b_k \xi(i-k), ~~~b_0=1.$$ (П68)

Здесь  $x(i)$  -- анализируемый  процесс c периодом дискретизации $T$, $\xi(i)$ -- белый шум с нулевым средним значением $\langle \xi\rangle =0$ и корреляционной функцией $C_{\xi \xi}(\tau)=\sigma ^{2}_{\xi } \delta(\tau)$ (см. Приложение). Спектральная плотность белого шума $S_{\xi \xi}(f)$ принимает постоянное значение для всех частот, равное дисперсии: $S_{\xi \xi}(f)=\sigma ^{2}_{\xi }$. Если все коэффициенты $b_k$ второго слагаемого равны нулю (кроме $b_0$, $b_0=1$), то получается так называемая модель авторегрессии (АР), т.е. модель процесса, значение которого в момент времени $i$ формируется как комбинация значений этого процесса в предшествующие моменты времени $(i-k),~~k=1,2 ...$ и значения независимой случайной составляющей $\xi(i)$. Если равны нулю все коэффициенты $a_k$, то модель называется моделью скользящего среднего (СС) и описывает безинерционное (не зависящее от предыдущих моментов времени) преобразование случайного процесса. Целые положительные числа $p$ и $q$ задают порядок моделей АР($p$) и CC($q$), соответственно. Когда $a_k$ и $b_k$ отличны от нуля, получается комбинированная модель -- авторегрессии - скользящего среднего АРСС($p,q$). Выбор типа модели предшествует спектральному анализу и делается на основе сведений о физической природе анализируемого случайного процесса. Кроме того, для выбора модели требуются предварительные сведения о возможной форме спектральной плотности мощности. Если исследуемый сигнал по предположению имеет спектр с острыми пиками, но без глубоких впадин (нулей), то наиболее подходящей является АР-модель, поскольку она дает спектр с резкими всплесками на отдельных частотах (рис. 15,  а). Если, наоборот, предполагаемый спектр имеет глубокие нули, но не характеризуется острыми пиками, то подойдет СС-модель (рис. 15,  б). АРСС-модель может применяться в обоих предельных случаях (рис. 15,  в). В тех случаях, когда одинаково пригодны несколько моделей, следует использовать ту из них, которая имеет наименьшее число параметров или, другими словами, наименьший порядок, поскольку оценки с хорошими статистическими свойствами можно получить, как правило, когда число оцениваемых параметров минимально. Кроме того, выбор модели зависит от имеющихся вычислительных ресурсов: АР-модели требуют значительно меньше вычислительных затрат, чем СС- и АРСС-модели. При неверном выборе модели оценка спектра будет ошибочной и, скорее всего, будет обладать худшими свойствами, чем оценка, полученная с помощью непараметрических методов.

:   Рис. 15. Примеры характерных функций спектральной плотности $S_{xx}(f)$ для АР-модели ( а), СС-модели ( б), АРСС-модели ( в).

Определим связь спектра последовательности $x(i)$ и параметров модели (1.68). Для этого применим преобразование Фурье (1.33) к выражению (1.68), используя свойство (1.45):

$$\begin{eqnarray*} F[x(i-k)]=X(f)e^{-j2\pi f k T}. \end{eqnarray*}$$

В результате получим следующее соотношение:

$$\begin{eqnarray*} X(f)=-X(f)\sum_{k=1}^{p}a_k e^{-j2\pi f k T} +\Xi(f)\sum_{k=0}^{q}b_k e^{-j2\pi f k T}, \end{eqnarray*}$$

где $\Xi(t)$ -- преобразование Фурье белого шума.

Перегруппируем  слагаемые  и выразим  спектральную функцию $X(f)$:

$$X(f) = \Xi(f)\frac{A(f)}{B(f)},$$ (П69)

$$A(f) = 1+\sum_{k=1}^{q}b_k e^{-j2\pi f k T},$$

$$B(f) = 1+\sum_{k=1}^{p}a_k e^{-j2\pi f k T}.$$

Поскольку для случайных процессов преобразование Фурье не определено, то необходимо перейти к спектральным плотностям. На основе соотношений (1.64) и (1.69) связь спектра мощности $S_{xx}(f)$ линейной стационарной модели (1.68) с параметрами модели и дисперсией белого шума запишется в виде:

$$S_{xx}(f)=\sigma ^2_{\xi} \left\vert\frac{A(f)}{B(f)}\right\vert^2.$$ (П70)

Таким образом, зная параметры $a_k$, $b_k$ и величину дисперсии шума $\sigma ^2_{\xi}$, можно определить спектральную функцию процесса.

Рассмотрим на примере авторегрессионной модели первого порядка ($p=1$) алгоритм нахождения спектральной функции анализируемого временного ряда с шагом дискретизации $T=1$. Модель (1.68) для $p=1$ имеет вид:

$$x(i)= - a_1 x(i-1) + \xi(i).$$ (П71)

В соответствии с выбранной моделью выражение (1.70) упрощается и принимает форму:

$$S_{xx}(f)=\frac{\sigma ^2_{\xi}}{1+2 a_1\cos 2\pi f+a_1^2}.$$ (П72)

Для того, чтобы оценить спектральную плотность $S_{xx}(f)$, необходимо предварительно найти оценки дисперсии шума $\sigma ^2_{\xi}$ и параметра $a_1$. Определим их. С этой целью сначала возьмем математическое ожидание от обеих частей уравнения (1.71):

$$E[x(i)]=-a_1E[x(i-1)]+E[\xi (i)].$$ (П73)

Математическое ожидание в первом слагаемом правой части (1.73) равно $E[x(i-1)]=E[x(i)]$, поскольку для стационарных случайных процессов математическое ожидание не зависит от начала отсчета времени. Второе слагаемое в правой части равно нулю вследствие того, что белый  шум  имеет нулевое среднее значение. Итак, выражение (1.73) имеет вид:

$$\begin{eqnarray*} E[x(i)]=-a_1E[x(i)]. \end{eqnarray*}$$

Это равенство  выполняется   только  при  следующем   условии: $E[x(i)]=0$. Это означает, что в модели авторегрессии первого порядка среднее значение процесса должно быть равно нулю.

Подсчитаем теперь дисперсию от обеих частей выражения (1.71):

$$\begin{eqnarray*} E[x^2(i)]= a_1^2 E[x^2(i-1)] + \sigma ^2_{\xi}. \end{eqnarray*}$$

Учитывая, что для стационарных СП $E[x^2(i)]=E[x^2(i-1)]$, получим соотношение

$$\sigma^2_x= \frac{\sigma ^2_{\xi}}{1-a_1^2}.$$ (П74)

Поскольку $\sigma^2_x$ -- величина положительная, на $a_1$ накладывается условие: $\vert a_1\vert<1$. Знание только оценки дисперсии не позволит оценить оба неизвестных параметра $a_1$ и $\sigma ^2_{\xi}$. Поэтому дополнительно подсчитаем автокорреляционную функцию $C(\tau )$ для времени задержки $\tau =-1$. В соответствии с формулой (1.31) для временного ряда $x(i)$ бесконечной длины ( $i=1, \dots , N$, $N\rightarrow \infty$) имеем

$$C_{xx}(-1)= \lim_{N\rightarrow \infty}\frac{1}{N}\sum_{i=2}^{N}x(i)x(i-1).$$ (П75)

Здесь учтено, что среднее значение процесса в модели авторегрессии равно нулю.

Далее подставим  в (1.75)   функцию  $x(i)$,   описываемую   формулой (1.71):

$$C_{xx}(-1) = \lim_{N\rightarrow \infty} \frac{1}{N} \sum_{i=2}^{N}(-a_1x(i-1)+\xi (i)) x(i-1) =$$ (П76)

$$= -a_1 \left[ \lim_{N\rightarrow \infty}\frac{1}{N}\sum_{i=2}^{N}x^2(i-1)\right] +\lim_{N\rightarrow \infty}\frac{1}{N}\sum_{i=2}^{N}\xi (i)x(i-1).$$

Множитель, заключенный в квадратные скобки, в соответствии с формулой (1.12) для бесконечного ряда ( $N\rightarrow \infty$) и нулевого среднего представляет собой дисперсию $\sigma^2_x$. Второе слагаемое правой части (1.76) равно нулю в силу статистической независимости процессов $x(i)$ и $\xi(i)$. Поэтому окончательно получаем выражение:

$$\begin{eqnarray*} C_{xx}(-1)= -a_1\sigma_x^2. \end{eqnarray*}$$

Из этого соотношения следует, что параметр $a_1$ равен:

$$\begin{eqnarray*} a_1=\frac {-C_{xx}(-1)}{\sigma_x^2}= -c_{xx}(-1), \end{eqnarray*}$$

где $c_{xx}(-1)$ -- значение коэффициента корреляции $c_{xx}(\tau)$ (см. формулу (1.25)) для $\tau =-1$.

Следовательно, искомый коэффициент АР-модели $a_1$ можно оценить по анализируемому процессу $x(i)$ конечной длины $N$ следующим образом:

$$\hat {a}_1=-\hat {c}_{xx}(-1)=\frac {-\sum_{i=2}^{N}x(i)x(i-1)}{\sum _{i=1}^N x^2(i)}.$$ (П77)

Для оценки дисперсии шума можно теперь воспользоваться выражением (1.74) и оценкой дисперсии (1.12):

$$\hat {\sigma} ^2_{\xi}=(1-\hat {c}^2_{xx}(-1))\frac {1}{N-1}\sum_{i=1}^N x^2(i).$$ (П78)

Вычисленные оценки параметра $a_1$ и дисперсии шума $\sigma ^2_{\xi}$ теперь можно использовать для оценки спектра процесса по формуле (1.72).

Как видно из рассмотренного выше примера, вычисление параметров авторегрессионной модели сводится к оценке моментных функций анализируемого процесса. Точно также с помощью моментных функций оцениваются параметры АР-моделей более высоких порядков. Оценка параметров СС-моделей является более сложной вычислительной процедурой, связанной с решением экстремальной задачи методом наименьших квадратов [7]. При вычислении параметров АРСС-модели задача еще более усложняется, поскольку сначала надо выбрать двумерный порядок модели $(p,q)$, затем найти параметры АР-части модели и, наконец, параметры, относящиеся к СС-части. Вычислительные сложности при нахождении параметров моделей связаны с необходимостью решения систем линейных (а в некоторых случаях -- нелинейных) алгебраических уравнений большой размерности. На сегодняшний день существует большое количество статистических пакетов прикладных программ, эффективно решающих задачу параметрического спектрального оценивания и помогающих оценить ошибку выбора порядка модели. Однако решение главного вопроса о выборе типа модели остается за исследователем.