Применение модели авторегрессии и скользящего среднего для спектрального оценивания в Фурье-спектроскопии Текст научной статьи по специальности «Математика»

В.А. Вагин, И.Б. Рейзин

ПРИМЕНЕНИЕ МОДЕЛИ АВТОРЕГРЕССИИ И СКОЛЬЗЯЩЕГО СРЕДНЕГО ДЛЯ СПЕКТРАЛЬНОГО ОЦЕНИВАНИЯ В

ФУРЬЕ-СПЕКТРОСКОПИИ

Традиционная процедура восстановления спектра (нахождения спектральной оцен-

• • • * * « 1

ки) В Фурье-спектроскопии опирается на теорему Винера-Хинчина [1], связывающую спектральную плотность мощности исследуемого процесса с его автокорреляционной функцией (в нашем случае - с интерферограммой) посредством преобразования Фурье, Соответственно для интерферограммы, дискретизированной согласно теореме отсчетов [2] , справедливо соотношение

со

Б (V) «НЕ 1<п)бхр<-12тгоп1|),

п=-<»

где - интересующая экспериментатора спектральная плотность мощности иссЛе-

I . . •

дуемого излучения;

Ил) - отсчеты регистрируемой интерферограммы с оптической разностью хода, равной пН;

И - шаг дискретизации интерферограммы; V - волновое число на интервале [О, 1/2Н].

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

5М а И Е 1(п) ехр>(-12т>пЮ,

П = -оо

где

1С п) =<

Кп), I п|^N О , IпI^

N - число зарегистрированных отсчетов интерферограммы;

I - предел изменения оптической разности хода в интерферометре в процессе указанной регистрации.

й

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

обусловленных боковыми лепестками аппаратной функции, интерферограмму аподиэи руют, то есть умножают на А(п) - некоторую весовую функцию» Соответственно

А 00 . . >ч

Б (V) = И 2 А (п) I < п) ехр(-1 2пл>пЮ .

П = -оо

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

До настоящего времени увеличения спектрального разрешения в Фурье-спектро-скопии добивались практически лишь путем увеличения пределов изменения оптической разности хода в интерферометре, что позволяло уменьшить ширину аппаратной функции прибора. Разработчики шли по пути усложнения приборов, что сопровождалось увеличением их габаритов и стоимости.

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

Одним из возможных подходов к решению этой задачи является использование модели авторегрессии и скользящего среднего (АРСС) исследуемого процесса, описываемой следующим рекуррентным соотношением [з] '

Р Ч

X < п) + Е акХ(п-к) = Е Ь^п-П,

к = 1 1=о

(1)

где

Х(п) - дискретный стохастический процесс;

Е(п) - "белый1' шум с нулевым средним и дисперсией г)2;

р, ц - порядки авторегрессивной (АР) и скользящего среднего (СО частей процесса;

ак/ " действительные константы.

Спектральная плотность мощности процесса (1) может быть записана [3] как

Б (у) = П2Ь

3

Е Ь [^ехр(-1 I)

1 = о

--р-

1 + 2 аь. ехр(-т 2тсуНк) к = 1

(2)

Таким образом, в рамках предлагаемой модели спектр может быть восстановлен по конечному числу параметров [т\2, а ^ ( к = 1, . •. р) , ЬСI = 0, .. .д)} ♦

Для обоснования выбора модели (1) для оценивания спектра в Фурье-спектроско пии заметим, что АРСС процесс является регулярным процессом, удовлетворяющим условию Пэли-Винера [4] 1/2Ь

/ I 1пБ (V) I с1\> <«>.

-1/2И

Известно [5], что регулярный процесс является реакцией физически осуществимого фильтра, на вход которого подается стандартная некоррелированная последовательность .

Фильтр называется физически осуществимым (каузальным)/ если его импульсная переходная характеристика ЬСО удовлетворяет условию

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

Определим автокорреляционную функцию 1(1) процесса (1). С учетом (3) она примет вид

-2 akI(l-k) ♦ S bkIcv(l-k), |l|:Sq,

1(1) =

k = i P

-Б akI(l-k), k = i

k = l

(4a)

Ul>q,

(46)

где I-y(l) - функция взаимной корреляции процесса XСn> и шума Е(п).

Подсистему (46) называют уравнениями Юла-Уокера. Параметры ак являются решениями этой подсистемы- Из (^б) видно, что неизвестные отсчеты автокорреляционной функции представляются в виде линейной комбинации р ее предыдущих от-

9 . »

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

Преобразуем выражение (2) к более удобному виду. Считая aQ«1 и учитывая, что а^О при к<0 и к>р и Ь^=0 при КО и l<q получим

Е a^expi-i2nvhk) к = о

Р Е т=

Amexp(-i 2nvhm),

-Р

где

т

Р

2 а ьа х. к=° к

-т'

< 5)

Е bLexp(-i2uvhl) l = o

Я

Е

т =

B^expC-i 2nvhm),

-q

где

Bm = Е ЬLЬL_т.

L =°

Пусть Bm=T|2hB|Jj, тогда

с учетом четности функций Ат и Вт получим

S (v) =

Е Bmcos(2nvhm) т = о_

Р

Е AmcosC2Tivhm) т=о

(6)

Коэффициенты Вт найдем, умножив обе части последнего выражения на знамена тель правой масти и выполнив преобразование Фурье:

В

m

Р

Е п =

A nI(m-n)

-Р

Как видим» при использовании выражения необходимости определять параметры Г|а и Ь^

(7)

(6) для восстановления спектра нет Точность же спектральной оценки за-

висит от точности определения параметров а^.

Реальный сигнал, снимаемый с Фурье-спектрометра, всегда зашумлен. Основным источником шумов в ИК диапазоне является тепловой шум приемника излучения. Этот

шум хорошо описывается распределением Гаусса с нулевым средним и некоторой дисперсией о2. Без ограничения общности будем считать шум белым (в случае его коррелированности можно применить декоррелирующее преобразование). Дисперсия шума может быть измерена до проведения основного эксперимента.

Определим отсчеты интерферограммы как

1<1>=1*(1>-еП>, <8>

где

I*(i) - идеальные значения автокорреляционной функции; С(i > - значения шума, удовлетворяющие условиям:

М{е<1)}=0, M{e(i)a}=t?a, M{eCi)e(j>}«0 при î*j. С учетом (8) система уравнений (46) примет вид

Р Р

Kl) + Е akI(l-k) = e(l) + Е ake(L-k)/ |ll>q. (9)

k=i k=1

Пусть I меняется от д + 1 до а+1:, где 1>р.

Введем обозначение:

Р

У<1) = еС I) + £ аке(1-к)

к = 1 к

и перепишем систему (9) в матричном виде где

Т=[1(р + 1>,...1<р+0]Т;

^ " > / • ■ • Р

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

Найдем корреляционную матрицу вектора V

С|к_1|=М{У(к)7(1)}=а3[б(к-1) ♦ а(к.и + Д а .а , ,+.] ,

к/ I = а^=0 при т>р,

где 6(к-1) - символ Кронекера. Введем обозначение

<5 Са) « [б( к -1) + а .. .. + £ а,а.. , . . . 1 ь I- х1

I к-1I ^ =1 } Iк-1|+]лк,1 = ч + 1,..+ 1

C-=aaGCa> .

Тогда корреляционная матрица вектора V запишется в виде

В работе [6] описано использование авторегрессионной модели для отыскания спектра электронно-циклотронного излучения плазмы.

Но подобная модель соответствует довольно узкому классу спектров- Известно [3] , что АР модель хорошо описывает пики, а СС модель - провалы в спектрах. Объединение этих моделей позволяет существенно расширить класс восстанавливаемых спектров, что продемонстрировано на рисунках 1-3.

На рис. 1 представлен модельный спектр (с характерными для спектроскопии Формами линий, а именно, гауссовой, лоренцевой и их суперпозицией), на основе

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

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

Рис. 1. Идеальная модель спектра

Рис. 2. АРСС оценка спектра. Порядок р АР части - 21, порядок ц СС части - 70. Число использованных точек интерферограммы - 92

Б (отн.ед.)

т

7 V (отн.ед.)

Рис, 3. Оценка спектра, полученная традиционным методом Фурье-анализа по 92 точкам интерферограммы

Запишем функцию правдоподобия для вектора V. Поскольку V распределен по нормальному закону, то с учетом (10) и (11) эта функция примет следующий вид:

Р(У/а, р)

1

(2тхо2Л/21б(а)|1/а

ехр {--— (7*Иа)Т6(а)~1 (Т+1а)} ,

2 о

где

1б(а)| - определитель матрицы б(а),

р - порядок АР части процесса.

Для нахождения параметров необходимо максимизировать по параметрам логарифмическую функцию правдоподобия:

1

1п Р(7/а, р) = -51пС2тш2) - ^1п|6(а) |--1~(Т+1а)Тб(аГ1 (Т+1а) .

2с2

Максимизация последнего выражения эквивалентна минимизации следующего функ-

(12)

ционала:

ф(У/а, р) = 1п|6(а)| + — (Т+1а)6(а)~1(Т+1а).

(13)

<т

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

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

но упрощается .

Пусть 1к(П - к-я реализация Тогда система (9) примет вид:

интерферограммы, а = - шум к-й реализации.

1о(1) ♦ Е ак1к(1-к) = е0С1)

к = 1

+ 2 акек(1-к), II1>я к = 1

(14)

ном случае

Найдем корреляционную матрицу вектора

описывается выражением:

р

У(1) = е0(I) + Е акек<1-к).

к = ч

V, каждая компонента которого в дан-

Поскольку М {EiCk)ej(1)}-о при (разные реализации шума статистически не-

зависимы), а также М {еi С k)е j ( I)} »0 при Ml получим

Р „

С.. ,, = аа6(к-1><1 + Е а?) ik-U j = 1 J

или

^ = + aTl)Etxt/

где Е - единичная матрица размерности t*t. Другими словами, V - гауссов "бе

»J ш

лый11 шум с нулевым средним и дисперсией о2 ( 1 + а а) .

Обратим внимание на то, что система (14) использует не всю имеющуюся у нас информацию. Из N реализаций интерферограммы можно составить N систем уравнений вида

Im(l) ♦ z aklj(t-k) = ej(l) + Б a^U-k), i ll>q, m*0,,..N-1, K>p, (15)

0 k=1 k=1

<<i-k>: e<->BOd Nu-k,i k=0,...,p, m=0,...,N-1 .

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

_ Р

Vm(О = е™(L) + Е akEk(l-k), m=0,...,N-1.

k = i

Зта функция будет иметь вил

N-1 _

Е 4Tm .+ I.ai/CT. >• 1яа)

ф(V/a, р) - - ![]}£ 1п2по^(1 + аТа) + —-7- . (16)

N 2 2аа(1 + а а)

Для нахождения оценки вектора параметров а нужно решить систему уравнений

Р> . 0/ п»1.....р.

В результате дифференцирования получим

N-1 . - М-1 ^М-1 т ^

И.1 М-1 гшТ.Х.1Та1ш Г Х.1.а

т = о т=о 1 + а а

л

Искомой оценкой будет такой вектор а, при котором выражение (16) достигает глобального максимума.

Неизвестным остался порядок р авторегрессионной части АРСС процесса. Его нельзя определить по методу наибольшего правдоподобия, поскольку изменение порядка меняет вид плотности распределения. Воспользуемся методом максимума энтро пии [7]. Как известно, этот метод позволяет построить функцию распределения, наилучшим образом удовлетворяющую экспериментальным данным.

Выражение (16) является естественной оценкой математического ожидания логарифмической функции наибольшего правдоподобия, которое задается выражением

я(е0, е) = <У/е0) 1п f(v/в)dv = м{1п ну/е)}, (17)

00

где

во = [а°,..,,ар , Р0]Т - истинный вектор параметров;

1 о

в » [а1,...,ар, р]Т - его оценка.

В работе [7] показано, что для произвольных плотностей распределения Их)

и д(х) выполняется соотношение + 00 +00

-/ ПхНп д (х) с1х>-/ Пх>1п f(x)dx. (18)

— ао — со

Причем эти выражения равны только при Их)=д(х)т Правая часть выражения (18) есть не что иное, как энтропия случайной величины х „ В обозначениях (17) послед нее выражение перепишется как

-|?(в0, 9)>-ме0, в0) = н(е0),

где Н(в0) - энтропия случайного вектора V.

Таким образом, выражение (16), взятое со знаком минус, является оценкой энтропии вектора V.

По определению [7] мерой энтропии случайного процесса называют предел

Н

1

х

^Н(*1/в..х2) при

В нашем случае оценкой меры энтропии будет величина

N-1

Hj(в) = i

S СТВ Ima)T(lm * .1яа>

In 2tio2(1 + aTI) + —

N t а2 С1 + aTa)

(19)

Известно [7], что мера энтропии "белого11 шума, распределенного по закону Гаусса с дисперсией а2, определяется выражением

НХ = I(ln 2псг2 ♦ 1). <20)

Этой величиной определяется мера энтропии случайного шума е, наложенного на автокорреляционную функцию. Другими словами, выражение (20) определяет меру энтропии спектральной оценки, полученной традиционным методом Фурье-анализа • Сравнивая выражения (19) и (20), видим, что увеличение энтропии - это та цена, которую приходится платить за увеличение разрешения, то есть АРСС модель может эффективно использоваться для спектрального оценивания при небольших шумах.

Для определения порядка р может быть использован информационный критерий Акайке (ИКА) [8], согласно которому выбирается такой порядок р, что минимально выражение

ИКА(р) = Н j(в) + p/t. (21)

Второе слагаемое в (20) отражает тот факт, что при увеличении порядка про-

/ч

цесса увеличивается дисперсия оценки После определения оценки вектора пара-

Л

метров а и порядка р АР части АРСС процесса могут быть найдены коэффициенты Ат (5), Вт (7) и оценка спектральной плотности мощности (6).

Реализация предложенного метода восстановления спектра будет представлена

в последующих публикациях.

Литература

К X у р г и н Я.И., Я к о в л е в В.П. Методы теории целых функций в радиофизике, теории связи и оптике. М.: Физматгиз, 1962.

2. Ш е н н о н К. Работы по теории информации и кибернетике. Под ред. Р.Л. Добрушина и О.Б. Лупанова (пер* с англ.). М.: ИЛ, 19 6 3 •

3. К э д з о у Дж.А. ТИИЗР. 1982, т. 70, N1 9. с. 256 .

4. Papoulis A. SIAM Review. 1985, v. 27, N 3.

5. к о р о л ю к 1 B.C. и др. Справочник по теории вероятностей и математической статистики. М.: Наука, 1985.

6. I w a m а N., I п о u е A., Tsukishima Т. et al. J. Appt. Phys., 1981, v. 52, N 9, p. 5466.

7. P a p о и L i s A. IEEE Trans. Acoust., Speech, Signal Processing. 1981, v. ASSP-29, N 6, p. 1176.

8. A k a i k e H. IEEE Transactions on Automatic control, 1974, v. AC-19, N 6, p. 716.