УДК 621.396.663
М. Е. Шевченко
Санкт-Петербургский государственный электротехнический
университет "ЛЭТИ"
Алгоритм совместного обнаружения и оценивания параметров источников радиоизлучения
Представлен алгоритм совместного обнаружения и оценивания параметров множества источников радиоизлучений (ИРИ) и излучаемых ими сигналов при трехэлементной антенной решетке. Алгоритм относится к классу ESPRIT-алгоритмов, основанных на формировании сигнального подпространства из наблюдаемых данных. Представленный алгоритм, в отличие от существующих, позволяет оценивать форму сигналов ИРИ и их параметры в условиях неопределенности сигнально-помеховой обстановки. Приведены результаты исследования алгоритма по модельным и по реальным сигналам.
Радиомониторинг, источники радиоизлучения, , азимут, угол места, обнаружение, оценивание, форма сигнала, антенная решетка, сигнальное подпространство
Основная задача панорамного радиомониторинга состоит в оценке частоты и в оценке угловых координат (азимута и угла места) множества источников радиоизлучения (ИРИ). Интерес представляет также определение типа сигнала.
В настоящей статье представлен алгоритм совместного обнаружения и оценивания частоты, азимута, угла места и формы сигналов от множества ИРИ, сигнальные составляющие которых присутствуют в наблюдаемых данных. Алгоритм разработан на основе ESPRIT-подхода [1], [2] к задаче синтеза алгоритмов определения направлений ИРИ и принципа инвариантности [3] для преодоления априорной неопределенности уровня аддитивного шума. ESPRIT-Подход (Estimation signal parameters via rotational invariance techniques) заключается в разделении сигнального и шумового подпространств на основе принятых данных, а также в формировании требуемых оценок на основе сигнального подпространства с учетом инвариантности антенной решетки к сдвигу.
Модели антенной решетки и сигнально-помеховой обстановки. Антенная решетка (АР) (рис. 1) состоит из трех всенаправленных антенн, расположенных в вершинах равностороннего треугольника. Длины сторон треугольника равны А, причем Afmax/c < 0.5, где fmax
- максимальная частота из диапазона принимаемых сигналов. Виртуальная линия, соединяющая антенны Y и X, ориентирована по направлению "север-юг" (антенна Y расположена севернее антенны X ).
Пеленг 9 k отсчитывается против часовой стрелки между направлением на ИРИ и
направлением на север. Азимут 9ak отсчиты- "èx
Рис. 1
9k
%
А
60°
V
Z
60°
А
60°
N
А
S
Известия вузов России. Радиоэлектроника. 2009. Вып. 1======================================
вается по часовой стрелке от направления на север. Поэтому = -0к .
Наблюдаемыми данными являются выборки х = (Хд, ..., хN — ), у = (уо, . ••, У^-1)
z = ( гд, ..., ZN—1) , полученные дискретизацией процессов, принятых элементами АР X, 7,
Z соответственно, с частотой Р . Отсчеты выборок записываются в виде
d
х = х(1/Р)= 2 Ьч (>/Р) ехр [/2/ (1/Р)] + £х (¡/Р);
к=1 d
Уг = У (¡/Р) =2 ЬкЧ (¡¡Р) ехр [/2/ (1/Р)] ехр (/ухук) + £у (¡¡Р);
к=1 d
Ъ = г (¡Р) = 2 ЬкЧ (¡Р) ехр [/2/ (г/Р)] ехр (/уХ2к ) + £2 (¡/Р),
к=1
где d - неизвестное число ИРИ; Ьк, /к - амплитуда и несущая частота сигнала к-го
ИРИ соответственно; Sk = {% (г/Р), ..., % [(N - О/Р ]| - выборка из комплексной огибаюТ
щейSk (t) сигнала к-го ИРИ, причем |[Sk (t)] Л = 1; уХук , ухгк - фазовые сдвиги между
0
сигналами к-го ИРИ в антеннах X, 7 и X, 2 соответственно, обусловленные разностью хода между элементами АР; £ х, £ у , £ г - выборки аддитивного шума в антеннах X, 7 и 2
соответственно.
При синтезе алгоритма предполагалось выполнение условия узкополосности принимаемых радиосигналов, поэтому 8к ^-^Ху(хг)] = % (t), где тху(хг) - задержка распространения сигнала между антеннами X и 7 (X и 2 ). Фазовый сдвиг между ¡-м и (г + п) -м отсчетами, принадлежащими принятому сигналу от к-го ИРИ в каждой из выборок, равен 2/п/Р. Фазовые сдвиги ухук , ух2к выражаются через пеленг 0 к и угол места Р к к-го
ИРИ как уххк =(2л/с)[/кА cos (60 -0 к)cos рк ] и ухук =(2л/с)(/кА 0к cos рк). Выб°рки
из аддитивного шума
^х =
£х I Р ], £х [ Р
^ =
£ у I , £
( N -1 у I Р
\г =
£г 1 о ], ., £г [ р>
присутствующие в наблюдаемых данных, полагаются взаимно независимыми с гаус-совской маргинальной плотностью распределения с неизвестной дисперсией.
На основе такой постановки задачи с применением обобщенного критерия наименьших квадратов, ESPRIT-подхода, методов линейной алгебры и принципа инвариантности разработан алгоритм совместного обнаружения и оценивания параметров ИРИ, который оценива-
======================================Известия вузов России. Радиоэлектроника. 2009. Вып. 1
ет число ИРИ d, формирует совместные оценки частоты , пеленга 0£, угла места Р £ и формы сигнала %.
Алгоритм совместного обнаружения и оценивания. Включает следующие этапы. 1. Формирование матрицы наблюдений
7 =
R =
^0
X Y 7
X =
х0 х1 . . xN-n-l Уо У1 . . ^-п-1
х2 . . ^-п ; Y = У1 У2 . . yN-n
хп-1 хп . . ^-1 _ _ Уп-1 Уп . . УN-1
гЫ - п-1
^2 ... 2Ы-г
2п-1 2п
-1
; п > d.
2. Вычисление корреляционной матрицы Ж = R (R*) ( - знак комплексного сопряжения; Т - знак транспонирования).
3. Разложение матрицы Ж = Е diag (Л)( Е *) по собственным числам Л = (Х0, ..., Xd-1, ..., ^зп-1), X) >Хd-1 > ^3п-1 и собственным векторам Е = ^ E2 ... Eзn-1}.
4. Оценка количества ИРИ d сравнением собственных чисел с порогом Са = С'а^, определяемым уровняем вероятности ложной тревоги а (С'а - параметр, зависящий от размера выборки N и параметра алгоритма п; а^ - состоятельная и эффективная оценка уровня шума).
5. Выделение из пространства собственных векторов Е = (, Е^) сигнального подпространства Е5 = (Е^ Е^ Е^ )Т образованного векторами, соответствующим d собственным числам, превысившим порог Са.
6. Получение оценок угловых параметров ИРИ. Для этого требуется:
6.1. Сформировать матрицу =(Е^ Е^ Е^ ) и вычислить корреляционную матрицу
ЕЕц ЕЕ12 ЕЕ13 ЕЕ21 ЕЕ22 ЕЕ23 _ ЕЕ31 ЕЕ32 ЕЕ33 _ из которой создать матрицы
EExy =
ЕЕ = () Ез,
_ ЕЕ11 ЕЕ12 и EExz = " ЕЕц ЕЕ13
_ ЕЕ21 ЕЕ22 _ _ ЕЕ31 ЕЕ33 _
6.2. Найти матрицы собственных векторов
ЕГху =
ЕГхУП ЕГхУ12 ЕГху 21 ЕГху 22
и EГxz =
ЕГхг1 1 EГxz
12
ЕГхг 21 ЕГхг 22.
матриц ЕЕХу и ЕЕхг. Вычислить ¥ху =-ЕГху1о ЕГ
-1
хУ
ху12 ху22 ' ™
¥„ = -ЕГ_ ЕГ
-1
хг12 хг 22
. Раз-
ложить ¥ху = Т(ху)ГхуТ ^(ху) и ¥хг = Т(хг)ГхгТ ^(хг) по собственным числам Гху = diag(Гху^, ..., Гху^), Гхг = diag(Гхг^,...,Гхг^) и векторам, составляющим матрицы Т(ху) и Т(хг).
6.3. Сопоставить собственные числа ¥ ху и ¥ хг изменением порядка следования собственных векторов в Т (хг) при сохранении порядка следования векторов в Т (ху). Для этого
вычислить корреляционную матрицу ТТ = [т* (ху)] Т(хг). Определить Mi - номер столбца
максимального элемента каждой I = 0, Л -1 строки матрицы ТТ. Привести порядок следования собственных векторов в Т (хг) к порядку их следования в Т (ху), соответственно изменив и порядок следования собственных чисел Г'хг = diag(Г'хг^, ..., Г'хг^) , Г'хг^ =Гxzм¡
так, что ¥хг = Т(ху)Г'хгТ_1(ху).
7. Получение оценок частоты. Для этого требуется:
7.1. Создать матрицы Ш (Езх ) , и2(Езх ), и1 (Е5у ), и2(Е5у ), и1 (), и2(),
где
и1(и ) =
и
0,0
и
0,Л-1
и
п-2,0
и
п-2,Л-1
и и 2 (и ) =
и
1,0
и
1, Л-1
и
п-1,0
... и
п-1,Л-1
7.2. Вычислить корреляционные матрицы
Т
ии(Ея) = {[Ц1(Ея)]*} и2(Е„); ии (Еах ) = {[и1(Е^ )]*}Ти2 (Ег ) .
7.3. Разложить матрицы
ии (Ех ) = Т (Ея) (Ея) Т-1 (Е„); ии (Е*у ) = Т (Е*у )Ф (Е*у)Т - (Е*у)
по собственным числам
Ф (Е„ ) = diag [Ф (Ея)р ..., Ф (Е„ )Л Ф (Ег ) = ^ [Ф (Ег ):,..., Ф ( Ег )Л ]
и векторам Т (Е3х ) , Т (Еэу ) , Т (Еы ) .
ш (Еу ) = {[и1(Е*у )]*}Ти 2 ();
ии (Е^у ) = Т (Е^у )Ф( Е*у)Т-1 (Е*у);
ф(Е^у ) = ^ [ф(Еу )l, Ф( Еяу )Л
======================================Известия вузов России. Радиоэлектроника. 2009. Вып. 1
7.4. Сопоставить (аналогично п. 6.3) полученные наборы собственных векторов T (Esx ), T (Esy ) и T (Esz ) собственным векторам T (xy) и изменить порядок следования
собственных векторов T (Esx ) ^ T ( xy ), T (Esy T ( xy ) , T (Esz ) ^ T (xy ) .
8. Получение матриц Ф'(Esx), Ф'(Esy ) , Ф'(Esz ) в соответствии с измененным порядком следования собственных чисел и вычисление оценки частоты:
fk = median {fu [Ф'(Esx )], fk [Ф'(Ey )], fk [Ф'(Esz )]}
по трем независимым оценкам частот fk [Ф' (Esx)] = (P2я) arg [Фk (Esx)] ;
fk [Ф'(Esz)] = (P2я) arg ^k (Esz)].
9. Вычисление оценок азимута:
arg (fx k )" 05arg (f'xy k )
fk [Ф'(Esy )] = (P2Я)arg[Фk (Esy)] ;
Qak = -arctg и угла места:
ßk = arccos
VÖ75arg (f'xy k )
arg fxzk, - 0.5arg (f^k)_ 2 + arg2
°.75 ( 2KAfk/c )
k = 1, d.
10. Вычисление оценки формы сигнала: £ = (... ) = ЯТЕ*Т (xy) .
Результатами действия алгоритма являются совместные оценки (/ , Эа^, Рк , % ),
к = 1, ё, частоты, азимута, угла места и формы сигнала как функции временных отсчетов.
Исследование алгоритма. Результаты работы алгоритма получены для модельной сигнально-помеховой обстановки и обработкой записей реальных сигналов при размере выборки N = 4096 и параметре алгоритма п = 600 . Шум моделировался гауссовским случай-
2
ным процессом с нулевым средним и дисперсией а . В табл. 1 приведены параметры пяти
2 2 12
ИРИ с относительными частотами /о = /к/Р и отношением "сигнал/шум" ^ = Ь>к /а ,
к = 1... 5. При заданных значениях отношения "сигнал/шум" вероятность обнаружения превышает уровень 0.9 при вероятности ложной тревоги 0.001. Суммарный амплитудный спектр
Таблица 1
ИРИ fo 0a, О ß, О г дБ Сигнал
1 15.5 25 45 20 ЛЧМ
2 50.4 15 45 17 АМ
3 75.5 30 45 17 Немодулированный
4 90.5 45 45 17 То же
5 110.6 210 45 17 - II -
Ху
4.2
2.1
Цау Уул/
V
30 60 90
Рис. 2
fc
4
2
3
1
5
0
Таблица 2 хА-
ИРИ /о §а, Р, . ° 195
1а 22.6 27.0 36.17
1б 28.5 27.2 48
1в 29.5 25.7 39.3 130 - 1б
2 50.04 21 58.3 1а
3 75.7 31.3 54.2 65 _ \
4 89.9 42.1 47.4
5 109 211.3 31.34
\
4
30 60
Рис. 3
90
/о
наблюдаемых данных от трех антенн показан на рис. 2. Номерами указаны спектры сигналов, соответствующие каждому ИРИ.
В табл. 2 приведены сформированные алгоритмом оценки частоты, азимута и угла места. На рис. 3 изображены спектры оценок формы сигнала, соответствующие каждому ИРИ.
Из табл. 2 и рис. 3 видно, что для ИРИ1, излучающего ЛЧМ-сигнал, алгоритм сформировал три близкие оценки частоты, азимута, угла места. Также получены три оценки формы сигнала, в результате сложения которых образован сигнал, квадратурная составляющая которого изображена на рис. 4.
На рис. 5 приведена оценка мгновенной частоты (кривая 1) суммарного сигнала после низкочастотной фильтрации и заложенный при моделировании ЛЧМ-сигнала закон изменения частоты (кривая 2). Видно, что восстановленный закон изменения частоты практически
соответствует заданному1.
() Результаты работы алгоритма по за-
писям реальных сигналов. На рис. 6 приведен суммарный амплитудный спектр фрагмента записи длительностью 25.6 мс. Запись сделана в полосе 100 кГц, центральная частота 7040 кГц, частота дискретиза-
0
- 7 - 14
0
/о ЛЧМ
0.05
0.04
0.03 0.02
1050
2100 Рис. 4
3150
X,
2500
0
1050
2100 Рис. 5
3150
7.3 7.6
Рис. 6
/, МГц
5
3
0
I
о
0
7
I
о
1 На рис. 4 и 5 по осям абсцисс отложено относительное время - номер отсчета в выборке оценки сигнала, умноженный на отношение N1 (N - п) за счет изменения размера выборки. 28
еа,
200 -
100
7.3
Рис. 7
40 -
20 -
Ч /
7.3
Рис. 8
Яе (5)
- 50
- 100
Яе (5)
- 20
- 40
7.6
.П
7.6
/, МГц
/, МГц
/, мс
t, мс
Яе (5)
Яе (5)
- 55 - 110
t, мс
t, мс
Рис. 9
1
2
7
Р
О
1
0
7
б
а
в
г
Известия вузов России. Радиоэлектроника. 2009. Вып. 1======================================
ции 160 кГц. Расстояние между антеннами 10 м. Выполнены предварительная калибровка и коррекция амплитудно-частотных и фазовых характеристик трактов.
По результатам обработки получены частотно-азимутальная панорама (рис. 7) и панорама "угол места - частота" (рис. 8). Для 19 обнаруженных ИРИ сформированы оценки формы сигнала. В качестве примера на рис. 9, a—г приведены оценки формы сигнала для обнаруженных ИРИ. Параметры сигнала, представленного на рис. 9, а: fo = 7 003 581 Гц, 9a = 209°,
Р = 0°. На рис. 9, б представлена оценка формы сигнала с параметрами fo = 7 038 862 Гц, 9a = 218°, Р = 0°. Рис. 9, в представляет оценку сигнала, имеющего параметры fo = 7 056 095 Гц, 9a = 139°, Р = 0°. На рис. 9, г представлена оценка формы сигнала с параметрами fO = 7 038 977 Гц, 9a = 154°, р = 0° 2.
Из рис. 7 и 8 видно, что алгоритм различил сигналы двух ИРИ 1 и 2, близкие по частоте, но различающиеся по угловым координатам. Различие по частоте составляет примерно 115 Гц.
Точность всех оценок, в особенности оценки формы сигнала, зависит от параметра n алгоритма. Чем больше n, тем точнее формируется оценка. Однако с увеличением n существенно возрастают вычислительные затраты алгоритма.
С помощью представленного алгоритма можно формировать совместные оценки частоты, угловых координат и излучаемых сигналов ИРИ с различными несущими частотами. Увеличение количества элементов АР при соответствующем алгоритме позволит различать и формировать оценки ИРИ с одинаковыми несущими частотами.
Библиографический список
1. Roy R., Kailath T. ESPRIT-Estimation of signal parameters via rotational invariance techniques // IEEE Trans. on ASSP. 1989. Vol. 37, № 7. P. 984—995.
2. Lemma A. N., Veen van der A. J. Analysis of joint angle-frequency estimation using ESPRIT // IEEE Trans. on sign. proc. 2003. Vol. 51, № 5. P. 1264—1283.
3. Богданович В. А., Вострецов А. Г. Теория устойчивого обнаружения, различения и оценивания сигналов. М.: Физмалит, 2003. 316 с.
M. E. Schevchenko
Saint-Petersburg state electroteknical university "LETI"
Algorithm of joint detection and estimation of a radio emission sources parameters
The ofjoint detection and estimation algorithm of set of radio emission sources parameters and their signals under three-element array is presented. The algorithm concerns to a class of ESPRIT-algorithms based on signal subspace formation from the observable data. Presented algorithm, unlike existing, allows to estimate the signals form of radio emission sources and their parameters in uncertainty signal/noise conditions. The simulation results and real signals processing by algorithm are included.
Radio monitoring, radio emission sources, azimuth, elevation, detection, estimation, signal form, array, signal subspace Статья поступила в редакцию 16 марта 2009 г.
2 Для удобства регистрации сигналы перенесены в низкочастотную область спектра.