Динамический алгоритм разрешения неоднозначности в фазовом угломере космической системы определения местоположения наземного источника радиоизлучения Текст научной статьи по специальности «Физика»

УДК 621.396

В.И. Тисленко, А.А. Савин

Динамический алгоритм разрешения неоднозначности в фазовом угломере космической системы определения местоположения наземного источника радиоизлучения

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

Постановка задачи

Проблема несанкционированного доступа в системы спутниковой связи с использованием искусственного спутника Земли (ИСЗ) на геостационарной орбите (ГСО) решается путем определения координат или пеленгации незарегистрированных земных станций (НЗС) [1-3]. При этом требуется высокая точность определения местоположения НЗС, которую может обеспечить фазовый метод. Необходимость измерения пространственного угла на НЗС требует применения фазированной антенной решетки (в общем случае объемной).

На рис. 1 показана геометрия взаимного расположения ИСЗ и НЗС в гринвичской системе координат (ГСК) [4]. Полагаем, что НЗС расположена на поверхности Земли и неподвижна относительно нее. ИСЗ расположен на идеальной ГСО, то есть неподвижен относительно НЗС.

Рис. 1 — Геометрия решения задачи определения местоположения

При описании поверхности Земли ограничимся сферой с известным экваториальным радиусом Щ = 6378136 м.

Допустим, приемная антенная система выполнена в виде объемной антенной решетки из М антенн, на основе которых образованы (JVf-1) независимых баз. Каждой паре антенн соответствует вектор базы

dm=Ál-Áj, т = 1,2,...,(М-1), (1)

где А^ — вектор, определяющий положение i-й антенны в орбитальной системе координат (ОСК) [4]. Полагаем, что система координат, в которой производятся измерения, совпадает с ОСК. Обозначим ф0 и А0 широту и долготу НЗС. Волновой вектор электромагнитного поля в месте приема в ОСК имеет вид

К (ф0. ¿-о) = rr [cos («)cos (Р)sin («) cos (Р)]Г , (2)

где а = 4 (<Ро До) и Р =(фо»^о) — азимут и угол места НЗС в ОСК. Отметим, что для идеальной ГСО и неподвижной НЗС эти величины постоянны во времени.

В общем случае разность фаз сигналов, принятых i-й и j-я антеннами, равна

Дфт (f) = фг (i) - ф, (t) = (£(Фо Д0),4 (0) + Пт (f) , (3)

где пт (t) — случайная ошибка измерения разности фаз; (•,•) — скалярное произведение векторов, указанных в скобке.

Зависимость dm (t) в выражении (3) указывает на возможность изменения во времени ориентации вектора т-й базы под действием управления ориентацией ИСЗ.

Далее рассматривается задача точной оценки координат НЗС ф0 и Х0 на основе использования фазовой информации в ситуации, когда образованная система баз и направленные свойства приемных антенн не обеспечивают формирование однозначных и точных оценок угловых координат НЗС. В качестве информативных сигналов, подлежащих обработке с целью определения оценок ф0 и А,0, зададим выходные сигналы (М-1) пар фазовых детекторов, подключенных к соответствующим парам антенн, образующих базы.

На рис. 2 схематично показан график функции Дф (t), определяемой выражением (3); сигнал на выходе измерителя разности фаз (фазометра) s(t); сигналы на выходе двух фазовых детекторов (ФД) cos [Дф (<)] и sin [Дф (í)J , подключенных к выходам приемно-усилитель-ных устройств (ПУУ) в каналах с i-й и /-й антеннами.

(П

ПУУ

Дф(/)

_ Т-

. Фазометр I__,

---г—' I

ПУУ

X

л/2'

ФД

созДф(г)

ФД

в вычислитель

ятДф(?)

fWr fW

Рис. 2 — Обобщенная структурная схема измерителя

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

Задача фильтрации состоит в формировании оценок гринвичских сферических координат НЗС по совокупности сигналов с выходов 2(М-1) фазовых детекторов.

Оптимальный алгоритм фильтрации координат

Синтез алгоритма оценки координат НЗС выполним на основе методов теории нелинейной марковской фильтрации [5, 6]. Для представления алгоритма обработки в форме, удобной для реализации в цифровом устройстве, введем дискретное время к = кТ, где /г = 1,2,... — номер отсчета; Т — интервал дискретизации.

Введем вектор состояния системы «ИСЗ - НЗС» х = [ф0 А.0]Г. По условию задачи он удовлетворяет разностному уравнению

х {к) = [ф0 (к)Х0 (к)]т = х(к -1). (4)

Начальные условия системы (4) определяет случайный вектор х (0). Полагаем, что для этого вектора можно задать априорную плотность распределения вероятностей

(ПРВ) 1У"[х(0)]» математическое ожидание £[х(0)] = тх и дисперсионную матрицу Е (¿(о)-т;с).(г(о)-*ут

= Vx. Совокупность сигналов, подлежащих обработке, объеди-

ним в вектор наблюдений

т=

~нс [г (*), я (*), *]"

г* (А) hs [х(ft),п(ft),ft]

(5)

где zc (ft) = {k)...zcm (k)...2cM^ (A)JT , ¿s(ft) = [2ls(fe)...Z;(A)...2lf_1(ft)]T, причем ^ (ft) =

= cos [(* [5 (ft)], 4 (ft)) + nm (ft)] и 4 (ft) = sin [(* [x (ft)], dm (ft)) + nm (ft)] - статистически зависимые сигналы на выходе I-й пары квадратурных фазовых детекторов; л (ft) = {/tm (ft)} —

(М-1)-мерный вектор случайных дискретных гауссовых последовательностей фазовых ошибок с нулевым средним и ковариационной матрицей

£[n(ft1)-/iT(A^)] = Vn-5(ftl-ft2),

(6)

где 5(») — дельта-функция Дирака; Уп — ковариационная матрица, учитывающая взаимную корреляцию фазовых ошибок между каналами в совпадающие моменты времени. Корреляция возникает вследствие влияния собственных шумов того приемного канала, который является общим при образовании измерительных баз.

Известно [5], что оптимальная байесова оценка состояния при квадратичной функции потерь является условным средним значением состояния х(к) по апостериорной ПРВ

w\x(k)\Z(k)\, то есть

Цк)= {...]" х{к)-Уг[х{к)\г(к)\йх, (7)

XI Хп

где 2 (к) = {г (г)}*^ — последовательность наблюдений до момента времени к.

В случае, когда апостериорная ПРВ имеет гауссов вид, оператор (7) строго приводится к явным рекуррентным уравнениям для оценки и дисперсионной матрицы ошибки фильтрации, которые определены алгоритмом работы линейного дискретного фильтра Калмана [5, 6].

Особенность поставленной задачи в том, что уравнения (5) существенно нелинейны по отношению к вектору ¿с(А) и нелинейно связаны с вектором возмущений й(А). Для решения задачи фильтрации применим численный метод Монте-Карло [7, 8] расчета апостериорной ПРВ на конечном множестве выборочных точек (частиц) с последующим вычислением оценки (7) в виде выборочного среднего.

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

Щн)]

£йГхЮ(*)1

/=1 L J

= |p)(ft)4^)(ft)]},

(8)

где ненормированный весовой коэффициент частицы с координатами дг ' (А) равен

Д [*(г> (А)] = Д [*«(к -1)] [г (А) | 5М (А)].

(9)

Начальное распределение (й = 1) частиц в пространстве состояний и их весовые коэффициенты определены априорной ПРВ

[г (н) | г (к -1)] = ту [х (1) 12 (о)] = ш[х (о)]. (Ю)

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

^[г(о)]-Хй[г(|,(о)]в[г(о)-х«(о)], ц[хМ(о)] = ± , = х,...,Р> (ц)

где х(0) = [фо Xof; *«(0) = ГФ[,£) Л«'

весовой

- координаты узлов сетки; ц коэффициент /-й точки; Р = Р1- Р2 — общее число точек, зависящее от количества разбиений

по широте Р1 и долготе Р2 ■

Апостериорная ПРВ истинного положения НЗС формируется в виде

мг[х(к)\г(к)]~ ¿ц[г«(Л)]6[5(А)_гО(А)]. (12)

t=i

При поступлении k-ro измерения апостериорная ПРВ рассчитывается в виде совокупности обновленных в соответствии с выражением (9) коэффициентов. Необходимая для этого функция правдоподобия W z (k) | х^ [к) J = L (fe)J с учетом гауссовой совместной ПРВ фазовых

ошибок, многозначности и статистической зависимости функций cos(") и sin(•) имеет следующий общий вид:

L [*(i) (Л)] = Ly [zc (*)] Í2 [гс (к), x(i) (к), ft] Ig [zs (к), x(i) (к), к

причем

м-1

■**[>'(»)> П I r l2'

Li zc{k),x^{k),k = £c-exp

r=l

(13)

(14)

(15)

где число ветвей многозначности cos (•) R = 2м 1; константа С =-—--

(>/Гя) " y]det\n

qir) _ qW . arccos^(/г) J ...arccosj^...arccos^2^_j (A) Jj ; матрица Q^ = diag[Sj (r) ...

вектор

... Sm(r)... SM_j (г)] с элементами

Sn (r) = sign i sin

2-я

о M-m

(r-8)

, 0 < 5 < 1;

с элементами

вектор Й«= [в?* (**'>(*),*) ... (*<'>(*),*) ... ^(¿«(йЦ

^ (*).*) - (*)] А (*)) ,

равными истинной расчетной разности фаз сигналов, приведенной к интервалу [-тг; я];

¿(0 — волновой вектор (2), соответствующий точке с координатами х^ (/г); V"1 -обратная ковариационная матрица фазовых ошибок; функция

Is[zs(ft),#(ft),ft]=;£ J ... J ... J N(y-J%Vn)dy, (16)

Г=1<Н(Г) am (r) aM-l(r)

где гауссова (M—1)-мерная плотность вероятностей

ы(у; й<'1У„) = С • ехрГ-| • [у - й^ • V;1 • [у - й^

Пределы интегрирования в уравнении (16) учитывают сочетание знаков координат вектора zs (ft) и определяются следующими соотношениями:

О» (г) = §{sign[4 (ft)] -1} + 7i[sm (г) - sign[4 (ft)]} , bm (r) = am(r) + n .

Таким образом, алгоритм оптимальной фильтрации полностью определен выражениями (8), (9) и (13)-(16). Текущая оценка вектора состояния на ft-м шаге вычисляется по выборке из Р точечных масс (частиц) на каждом временном шаге в виде весового среднего (8). Причем при ее вычислении учитывается предыдущая информация (в виде весов и оценок на (ft-l)-M шаге) и вновь поступившие наблюдения z (ft) в виде правдоподобия каждой частицы (13).

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

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

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

Динамическое устранение неоднозначности выполняется за счет управления положением антенн ИСЗ, которое предполагает вращение ИСЗ вокруг оси, направленной на центр Земли.

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

Условия моделирования и результаты расчетов

Математическое моделирование алгоритмов выполнялось в программе, разработанной в среде Matlab 7, на ЭВМ с характеристиками: CPU Athlon 64 3000+, RAM 2 х 256 Мб Dual Channel. Исследование ошибок определения координат выполнено методом статистических испытаний при усреднении по ансамблю из 20 независимых реализаций.

Антенная решетка фазового пеленгатора состоит из четырех антенн. Первая база образована антеннами А1гА2 и параллельна экваториальной плоскости. Вторая образована антеннами А3 и А4 и параллельна меридиональной плоскости [1] (рис. 3). Третья база образована антеннами Аг и А3. Координаты антенн в орбитальной (измерительной) системе координат в длинах волн К заданы следующими:

А^-ЮОЛ; 500А.; 0), Aj(-100X; - 500Л; 0), А3 (0; 0; -250Я), А4 (0; 0; 250Л).

Отметим, что вектора баз не лежат в одной плоскости, то есть решетка является объемной.

Размер априорной зоны расположения НЗС на первой итерации составляет ±3° по широте и долготе относительно подспутниковой точки (точки стояния), имеющей нулевые широту и долготу. В действительности размер зоны ограничен диаграммами антенн. На расстоянии около 36000 км эта область на поверхности Земли ориентировочно имеет диаметр 600 км. Положение НЗС в этой зоне полагается случайным с равномерным законом распределения вероятностей. Расчет выполнен для трех значений СКО фазовой ошибки оДф равных 10, 20 и 30 градусов. Общее число точек Р = 10000 (то есть Рг = Р2= №0). Расчет выполнен для трех значений скорости вращения: 0,5; 1; 2 градуса за время, равное периоду поступления данных Т. На рис. 4 показана апостериорная ПРВ для четырех моментов времени.

Рис. 4 — Эволюция апостериорной ПРВ (Од<р = 20° , скорость вращения ИСЗ 1 °/Т ): а — г = Т; б — г = 2Т; в — * = 4Т; г — г = 8!Г

На второй итерации размер зоны составлял ±20'. Результаты моделирования приведены в таблице (через знак «/ » обозначены результаты после первой и второй итераций). Отметим, что при отсутствии относительного перемещения ИСЗ и НЗС точность оценок координат зависит не от общего времени измерения, а от количества N некоррелированных во времени отсчетов наблюдений.

Параметр Скорость вращения ИСЗ, °/Т СКО ошибки по разности фаз, градус

10 20 30

СКО оценки ф0, угл. мин. 0,5 1,34 / 1,1 2,2 / 1,69 2,75 / 3,98

1 1,01 / 0,93 1,57 / 1,44 1,42 / 2,07

2 0,99 / 0,68 0,94 / 0,6 1,24 / 1

СКО оценки А.0, угл. мин. 0,5 0,91 / 0,5 1,22 / 0,6 19,7 / 19,8

1 1,11 / 0,67 1,03 / 0,52 1,08 / 0,98

2 1,19 / 0,48 1,12 / 0,59 1,13 / 0,57

СКО места, км 0,5 2,94 / 2,1 3,6 / 3,27 36,7 / 37,3

1 2,72 / 1,9 3,41 / 2,79 3,25 / 4,19

2 2,81 / 1,51 2,66 / 1,53 3,06 / 2,1

Оценка времени, ед. Т 0,5 25 30 >30

1 15 20 25

2 10 10 15

Заключение

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

Наличие суточных эволюций ИСЗ на реальной ГСО также позволяет решать задачу устранения многозначности измерений, причем без вращения ИСЗ. При этом требуется значительно большее время. Использование ИСЗ на высокоэллиптической орбите позволяет сократить это время. В предложенном методе при заданной структуре антенной системы и СКО фазовых ошибок (10-30) градусов скорость сходимости оценок определяется временем, необходимым для поворота ИСЗ на угол не менее 12-30. При этом требуется 10-35 независимых измерений.

Литература

1. Панько С.П. Фазовая пеленгация в спутниковой связи / С.П. Панько, В.В. Сухотин / Электронный журнал «Исследовано в России». — С. 380—388. — Режим доступа г http:// zhurnal.ape.relarn.ru/articles/2003/035.pdf

2. Панько С.П. Несанкционированный доступ в системы спутниковых коммуникаций / С.П. Панько, В.В. Сухотин / Успехи современной радиоэлектроники. - 2002. - № 4.

3. Могучев В.И. Доплеровская пеленгация земных станций через геостационарный спутник связи / В.И. Могучев / «Электросвязь». - 2003. - № 1. - ISSN 0013-5771.

4. Жданюк Б.Ф. Основы статистической обработки траекторных измерений / Б.Ф. Жда-нюк. - М. : Сов. радио, 1978. - 384 с.

5. Сейдж Э. Теория оценивания и ее применение в связи и управлении / Э. Сейдж, Дж. Меле ; пер. с англ. ; под ред. проф. Б.Р. Левина. - М. : Связь, 1976. - 496 с.

6. Ярлыков М.С. Статистическая теория радионавигации / М.С. Ярлыков. - М. : Радио и связь, 1985.

7. Doucet A. On Sequential Simulation-Based Methods for Bayesian Filtering. Technical report CUED / A. Doucet / F-INFENG / TR 310, Department of Engineering, Cambridge University, 1998.

8. Спанье Дж. Метод Монте-Карло и задачи переноса нейтронов / Дж. Спанье, Э. Гел-бард ; пер. с англ. ; под ред. А.Д. Франк-Каменецкого. - М. : Атомиздат, 1972.

9. Jeroen D. Hol. Resampling in particle filters / D. Hoi. Jeroen / Division of Automatic Control, Department of Electrical Engineering. Lin. univ. 2004.

Тисленко Владимир Ильич

Канд. техн. наук, доцент каф. радиотехнических систем ТУСУРа

Телефон: (3822) 41 38 92

Эл. почта: wolar@post.tomica.ru

Савин Александр Александрович Аспирант каф. радиотехнических систем ТУСУРа Телефон: (3822) 41 38 92 Эл. почта: saasavin@mail.ru

Tislenko V.I., Savin А.А.

A dynamic algorithm of the ambiguity resolution in a phase angle meter of space-based system intended for location of a ground radiation source

An estimation problem of geographical coordinates of unregistered ground stations from geostationary satellite has been solved for the conditions when the satellite array does not allow to get the unambiguous and exact readout of angular coordinates by signals received at one point time. A processing algorithm and research results of probability characteristics of coordinates estimation are presented.