Вариационный непараметрический метод оценки двухмерного пространственного спектра сигналов, отраженных от ионосферы Текст научной статьи по специальности «Электротехника, электронная техника, информационные технологии»

УДК 550.388.2

ВАРИАЦИОННЫЙ НЕПАРАМЕТРИЧЕСКИЙ МЕТОД ОЦЕНКИ ДВУХМЕРНОГО ПРОСТРАНСТВЕННОГО СПЕКТРА СИГНАЛОВ, ОТРАЖЕННЫХ ОТ ИОНОСФЕРЫ

© 2010 г. Н.Н. Шкрылев

Южный федеральный университет Southern Federal University

ул. Зорге, 5, г. Ростов-на-Дону, 344090, Zorge St., 5, Rostov-on-Don, 344090,

physdekan@sfedu.ru physdekan@sfedu.ru

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

Ключевые слова: разделение лучей, непараметрический метод, двухмерная спектральная плотность мощности, частотно-пространственно-временная корреляционная матрица.

The possibility of definition of high spectral resolution angular spectrum ofHF signals reflected from the ionosphere is shown through the use of signal processing variational technique.

Keywords: signal reconstruction, non parametric technique, two dimensional space spectrum, frequency-space-time correlation matrix.

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

Известен ряд общих приемов преодоления этой проблемы, среди которых можно отметить следующие:

1. Аналитическое разделение лучей, требующее многоканальной приемной аппаратуры с выровненными амплитудно-фазовыми характеристиками и сохраняющее работоспособность только при хорошем отношении сигнал/шум [1, 2].

2. Разделение многолучевого поля на отдельные составляющие лучи по доплеровскому сдвигу частоты [3]. Недостатком такого подхода является ряд необходимых условий: жесткие требования к стабильности частоты передатчика, необходимость длительного периода наблюдений за сигналом.

3. Подавление влияния многолучевости статическим усреднением измеренных на антенных элементах решетки фаз напряжений относительно опорного элемента [4, 5]. Как следствие, способ оказывается применимым только к длительным выборкам сигналов.

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

Метод оценки двухмерного углового спектра с высоким разрешением

Будем исходить из интегрального представления Крамера для пространственного спектра стационарного пространственного поля сигнала:

X(r) = —(e'krdZ(k) • (1)

(2л)3 i ' '

Здесь X (r) - стационарное пространственное поле; Z(k) - пространственный спектр стационарного поля; r - радиус-вектор пространственной точки; k -волновой вектор. Пространственный спектр стационарного пространственного поля - функция с некоррелированными приращениями, обладающая свойствами:

(Z(Ä) - 0), если(Х(») = 0; (2)

(Z(k)Z *(k')) = (|Z(k)2|)s(k-V) •

Тогда имеет место следующее спектральное представление пространственной корреляционной функции:

B(r-r') = (X(r) — X *(r ')) = —Ц- le-*(r-r,)dF(k), (3)

(2л-) •

где F (k) - пространственный спектр пространственной корреляционной функции B(r -г'). Причем, если пространственное поле X(r) не содержит дискретных пространственных компонент, то из (1)-(3) следуют соотношения, устанавливающие связь пространственного спектра стационарного поля и спектра пространственной корреляционной функции:

(Z(Ä)Z *(*')) = (2K)3ö(k-k')dF(k)d3k ; (4)

(|Z(k)\2^ = (2Ж)3dF(k).

Если пространственное поле X(r) содержит дискретные пространственные составляющие, т.е. компоненты вида Aelk°r, то оно может быть представлено в виде X(r) = Aelk°r + х(г), где х(г) - стационарный случайный процесс, причем (x(f)) = 0, а

А = ае"р - комплексная амплитуда дискретной пространственной составляющей. В свою очередь, про-

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

вид В(г,г) = а2вЛо(р~Г) + (х(г)х*(г)) .

В этом случае соотношения (2) изменяются и приобретают вид

{г(к)г * (к')) = (2л-)6 а28(к0 - к)8(к - к')<Лък<Лък' +

+ (2л:)3 8 (к - k)dF(k)d3k ;

= (2ят)3 a 2S(k0- k)d 3k + (2ят)3 dF(k).

(5)

Итак, пространственный спектр пространственного стационарного поля - функция с некоррелированными приращениями, связь между ним и спектром пространственной корреляционной функцией устанавливается соотношениями (4), если пространственное поле не содержит пространственных дискретных компонент, и соотношениями (5), если таковые присутствуют.

Перейдем от непрерывных представлений к дискретным пространственным представлениям. Будем считать, что датчики поля размещены в дискретных пространственных точках с радиус-векторами гп, где п = 1, а N - количество антенных элементов.

Пусть поле в зоне приема формируется в результате интерференции ./ лучей с волновыми векторами 1<г1 = 0.../ -1. Тогда скалярное пространственное

поле Х(г) в пространственной точке с радиус-векторам гп может быть записано в виде

J-1

= 14

(6)

Хи=Х(ги) = ^еКк^

7=0 7=0

где А7 - комплексный случайный стационарный процесс. Очевидно, что дискретное представление скалярного дискретного поля Хп является аналогом непрерывного представления (1) непрерывного скалярного пространственного поля Х(г). Причем А7 имеет смысл дискретного спектра дискретного пространственного стационарного процесса Хп. Очевиден и смысл 2(к) при предельном переходе от (6) к (1), а именно й2(к) имеет смысл суммы всех комплексных амплитуд пространственных гармоник А7, попадающих в бесконечно малый объем d к. Если не содержит дискретных составляющих, то А7 попарно не

2

коррелированы: ^А7= \а^>71=^871+ ^871, где

а2 - мощность пространственного шума поступающего с 7 -го направления. Если Х(г) содержит, например, дискретную пространственную составляющую с амплитудой а0 , то имеет место соотношение: 2

{А:А*^\А:\ Ч/=аК/+СГ.

Дискретное представление (6) может быть переписано в матричном виде Л' Е.I. где Хт = -^"о ■ ■ - ^1 -вектор-строка распределения пространственного поля

по апертуре антенной решетки; Ат = Ио- - '1 ■■■■',/-1 ~~ вектор-строка дискретного пространственного спектра, причем каждый дискретный индекс 7 соответствует определенному 7 -му волновому вектору к7, т.е. отвечает определенному направлению в пространстве. Прямоугольная Лгх./ матрица Е имеет вид

E =

eih0r0 'Vi

eikiro

'Vi

ik

J- iro

ik

eik0rN- i e'kirN - i

ik r

ik

J-irN - i

(7)

Теперь может быть построена пространственная корреляционная матрица как среднее статистическое внешнего произведения вектора X на эрмитово сопря-

тт

женный вектор X . Ее элементы с учетом некоррелированности спектральных амплитуд А7 имеют вид

'пт у = 01 = 0 х '

J-IJ-1

2

^ ^ EnjEmj\Aj I

|2 J-1 * .......... 4'/ = 2 Ет-Ет;\А;

7 =0/ = 0

Элементы пространственной корреляционной матрицы, как показано в [6], оцениваются одновременно с частотно-временной локализацией обнаруженного сигнала, а в качестве оценок для (ХХН^ можно использовать частотно-пространственно-временные матрицы 5о, найденные усреднением в частотно-временной области существования сигнала. Тогда, зафиксировав столбец с номером т матрицы 5о , получим матричное уравнение = Е11 .

VI т г и

Это уравнение поэлементно имеет вид

^О— = Е ¡Е ;\А;

j=о

(8)

(9)

Неизвестный вектор-столбец ит состоит из эле-

ментов Б,

m]

AJ

причем |f/m| = \А\2. Как следствие,

задача оценки пространственной спектральной плотности мощности (СПМ) сводится к решению системы уравнений (8) или (9). Модуль найденного комплексного вектора ит дает пространственную СПМ.

Заметим, что система (8) является недоопределен-ной и имеет N уравнений (по числу антенных элементов) и У неизвестных компонент вектора ит. Естественно, на практике при оценке пространственной СПМ сигналов, отраженных от ионосферы, всегда будет выполняется неравенство ./ > Л'. Методы решения уравнения (8) обсудим ниже, а сейчас остановимся на возможности повышения значимости получаемых оценок пространственных СПМ.

Частотно-пространственно-временная мат-

рица 5о содержит N независимых столбцов (по числу антенных элементов в решетке). Как следствие, изменяя в уравнении (8) индекс т от 0 до Ы—1, найдем N независимых решений и оценок простран-

e

e

e

2

ственной СПМ. Тогда значимая оценка

A

про-

странственной СПМ определяется по формуле: |~|2 1 N-1

И .

I 1 Н т=0

Оценка пространственной спектральной плотности мощности

В общем случае при 3 > N уравнение (8) имеет бесконечное множество решений. С помощью псев-

и м м \ дообратной Еп = А" (ААП )-1 матрицы может быть

найдено единственное решение с наименьшей нормой. Однако, как показали предварительные оценки в ходе компьютерного моделирования задачи, получаемое таким образом решение является статистически неустойчивым даже при наличии малых погрешностей в оценках частотно-пространственно-временной матрицы ^о . В свою очередь такие погрешности всегда имеются и обусловлены не только конечной точностью измерений, но и шумами естественного происхождения. Поэтому для решения поставленной задачи (8) используем вариационный метод, идеи которого изложены в работе [7]. Реализованный нами алгоритм является его обобщением на случай оценки двухмерного углового спектра с последующей оптимизацией способа численных оценок. Обнаружение сигнала и его пространственная локализация осуществлялись на основе алгоритма, основные моменты которого состоят в следующем. Модель принимаемого сигнала представляется в виде Л' = ЕII+IV, где 5" - измеренный вектор-столбец амплитудно-фазового распределения принимаемого сигнала по апертуре антенной решетки с размерностью, равной числу антенн N (столбец частотно-пространственно-временной матрицы 5о); W - вектор шума размерностью N; V

- комплексный вектор размерностью N пространственной СПМ принимаемого сигнала КВ источника; Е

- матрица оператора Фурье размерностью NxJ. Для псевдообращения записанного уравнения мы будем исходить из условия оптимизации:

II = vamJ0(U), где ./()(//) имеет следующий вид:

и

J-1

9 j -1 . ,2 —

J(U) = |S-EU2+ Z (|U7| +s)2 , где p<2 и Я -

J= О

скалярные параметры;

- lp -норма.

Опишем численное решение оптимизации поставленной задачи. В общем случае для решения этой проблемы мы использовали аппроксимацию 1р -нормы

при р<\: 1(и) = р--Еи\\22+'1^(\и]\2+£)2 , где е -

7=0

малая константа; V у - у -й элемент вектора V ; 3 -длина вектора V . 3(V) минимизируется методом

Квази-Ньютона: и(п+1) = V" +]Н~\и(-п))Ш(и(-п)), где п - номер итерации; у регулирует размер шага на

каждой итерации; V./ (//) - градиент весовой функции; Н - аппроксимация Гессиана, представимая в виде

HQJ) = ЕНЕ + Я diagl

p/2

(U,| + e)l~ P/2

Здесь diag() - это диагональная матрица, состоящая из диагональных элементов, представленных в скобках. После подстановки последнего выражения в предыдущее получим следующий итерационный алгоритм: [[(!:""' = 2АНБ, где использовано значение у = \. Алгоритм работает до тех пор, пока не

V (п+V (п)

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

■ < S, где

D(U(n)) -

U(п)

S - малая константа.

Для повышения быстродействия была применена

тт _I _I

лемма об обращении матрицы: (D + En Е)"1 = £>_1 +

+ {D~lEH)(ED~lEH +/)"1 +(Е£>)-1, где диагональная матрица MiagQ размерностью ./ х./ ; Е - известная матрица размерностью NxJ ; I - единичная матрица размерностью N х N.

Применительно к нашему случаю получаем уравнение в виде U(n+V> ^[Lr1 ~(D~lEH)(W~lEH +Е)-1 (ED_1)]x

х (EHS).

После преобразования окончательно получим следующее итерационное уравнение:

u(n+\) =D~1EH[I-(ED~1EH +iy1(ED~1EH )]S .

Представленный алгоритм был обобщен применительно к двухмерному пеленгованию многолучевых KB сигналов, отраженных от ионосферы. Оценочное время расчета при р = ОД для диапазона сканирования по пеленгу а от 0° до 36СР и углу места А от 10° до 70° на процессоре Intel с частотой процессора 1,6 ГГц составляет примерно 1,5 с.

Алгоритм использован в задаче двухмерного пеленгования источников декаметровых волн применительно к 16-канальному пеленгатору. Антенная решетка пеленгатора располагалась на поверхности Земли на площадке 100x100i , использовался 16-канальный приемник с выровненными амплитудно-частотной и фазо-частотной характеристиками; сигнал по промежуточной частоте каждого канала оцифровывался 16-канальным АЦП с частотой дискретизации 2,5 МГц.

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

Полученные результаты двухмерного пеленгования представлены на рис. 1-4.

б

Рис. 1. Результаты двухмерного пеленгования в условиях многомодового распространения на трассе Гамбург-Ростов (2193 км), частота 10100.8 кГц. Разделены нижние и верхние лучи, отраженные от F-областей ионосферы. Показаны изменения от времени азимута (а) и угла места (б)

На всех рисунках изображены зависимости азимута а (а) и угла места А (б) от времени t. Причем рис. 1 и 2 - трехмерные (по оси 02 отложены значения пространственной СПМ), рис. 3 и 4 - двухмерные.

На рис. 2а, б показаны результаты для трассы Измир-Ростов на частоте 10112 кГц (азимут 230°, дальность 1390 км). Разработанный алгоритм оценки углового спектра позволяет разделить лучи, отраженные от Е- и F-областей ионосферы. Е-отражениям соответствуют углы места в окрестности 14°, в свою очередь Г-отражениям - углы места 25°.

На рис. 1 а, б показаны результаты для трассы Гамбург-Ростов на частоте 10100,8 кГц (азимут 297°, дальность 2193 км). Разработанный метод позволяет разделить по азимуту и углу места парциальные лучи, соответствующие нижнему и верхним лучам, женным от F-области ионосферы.

На рис. 3 а, б показаны результаты для трассы Анкара-Ростов на частоте 8052 кГц (азимут 216,6°, дальность 982 км).

а

б

Рис. 2. Результаты двухмерного пеленгования в условиях многомодового распространения на трассе Измир-Ростов (1390 км, 230°), частота 10112 кГц. Разделены лучи, отраженные от Е- и F-областей ионосферы

На рис. 4 а, б показаны результаты для трассы Москва-Ростов (источник сигнала - станция точного времени РВМ) на частоте 9996 кГц (азимут 354,4°, дальность 960 км).

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

б

Рис. 3. Зависимость азимута (а) и угла места (б) от времени на трассе Анкара-Ростов (азимут 216,6°, дальность 982 км, частота 8052 кГц). Большие кружки - оценки классическим методом, точки - оценки, полученные по нескольким максимумам разработанным вариационным методом.

Проведенные исследования позволяют сделать выводы:

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

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

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

t, c

t ,c

б

Рис. 4. Зависимость азимута (а) и угла места (б) от времени на трассе Москва (РВМ)-Ростов (азимут 353,4°, дальность 960 км, частота 9996 кГц)

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

Выражаю особую благодарность д.ф.-м.н., с.н.с., проф. каф. радиофизики ЮФУ Г.Г. Вертоградову за неоценимую помощь в написании данной статьи.

Литература

1. Gething P.J.D. Radio direction finding and superresolution.

London, 1990. 365 p.

2. Марпл-мл. С.Л. Цифровой спектральный анализ и его

приложения. М., 1990. 584 с.

3. Афpаймович Э.Л. Интерференционные методы радио-

зондирования ионосферы. М., 1982. 198 с.

4. Вертоградов Г.Г., Кондаков Е.В. Уменьшение влияния

многолучевости на точность определения углов прихода интерферометрическим методом // Радиоконтроль: науч.-техн. сб. Ростов н/Д, 2000. Вып. 3. С. 31-38.

5. Вертоградов Г.Г., Кондаков Е.В. Уменьшение влияния

многолучевости на точность определения углов прихода интерферометрическим методом // Радиотехника. 2003. № 1. С. 86-90.

6. Пат. № 2150122. Способ определения двумерного пе-

ленга и частоты источников радиоизлучения.

7. Cetin M., Malioutov D.M., Willsky A.S. A variational tech-

nique for source localization based on a sparse signal reconstruction perspective // Processings of the 2002 IEEE International Conference of Acoustics, Speech and Signal Processing. Orlando, USA. 2002.

Поступила в редакцию

3 июля 2009 г.

а

а