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

УДК 621.396

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

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

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

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

Геометрия задачи определения местоположения показана на рис. 1, где ОХУ — декартова геоцентрическая система координат (ГСК). Ее центр расположен в центре Земли, ось ОХ находится в плоскости экватора и направлена в точку весеннего равноденствия; ось 02 параллельна оси вращения Земли и ориентирована в направлении севера; направление оси ОУ соответствует правой системе координат [1].

Рис. 1 — Геоцентрическая система координат

Текущее положение на орбите ¿-го КА (наблюдателя) в момент времени t определено заданием радиус-вектора roi (t), положение ИИ определяет радиус-вектор rg (i). Таким образом, вектора

Л1; (О = ro¡ (*)"г.(*) и Av¿ (t) = roi{t)~r3{t) = voi(t)-vs (t) (1)

определяют в момент времени t соответственно текущее взаимное положение и скорость направленной линии визирования «источник — наблюдатель» в ГСК.

Пусть ИИ имеет частоту радиосигнала f3 и расположен в точке на поверхности Земли, широта и долгота которой равны соответственно <ра и Xs. В гринвичской системе координат [1], вращающейся вместе с Землей, широту и долготу точки определяет постоянный во времени вектор гринвичских координат источника

Tsg = [^Е ' cos Ф8 "COSК&Е ■ C0SФа ' s*nК^Е " S"1 Ф8Г ' (2)

где Re = 6378136 м — радиус Земли. Таким образом, в гринвичской системе координат вектор скорости ИИ rsg = v$g = 0. Вектора положения и скорости в ГСК и гринвичской системе координат связаны известными [1] соотношениями

(О = [у. (О (0 z* (0]Г = и (0 • Гв, ; Ve (t) = Q£ • [-„, (t) *8 (t) 0]r, (3)

где матрица Х1(<) учитывает суточное вращение Земли; &Е= 0,729211-1(Г4 с-1 — средняя

угловая скорость вращения Земли.

Введем трехмерный вектор переменных состояния

=[*1 *2 *з] = [/. К Ф.], (4)

координаты которого постоянны во времени и являются неизвестными случайными величинами. Очевидно, что для неподвижного ИИ в дискретном времени для вектора х справедлива система разностных уравнений

х(А) = х(А-1), И = 1,2,.... (5)

Случайные начальные условия для уравнения (5) определены заданием априорной плотности распределения вероятностей (ПРВ) Ж(х(0)). Априорное математическое ожидание

= V,

о, имеет диагональ-

М [х(0)] = ш0. Ковариационная матрица М (х(0) - т0)(х(0) - т0)Т

(2 2 2 \

о01 о02 ст031. Далее по тексту М — оператор статистического усреднения.

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

х<м (г) = [гог (О уог • Известно, что в этом случае хы (<) является решением системы нелинейных дифференциальных уравнений. В векторной форме они имеют вид [1]

(0 .«О (*)_

(6)

где ас(г) = -ц-гог(г)/|го;(г)|3 — гравитационное ускорение и ц = 0,3986013 Ю15 кг м/с2 — произведение массы Земли на гравитационную постоянную. Начальные условия хо£ (0) для

системы (6) определяют орбиту 1-го наблюдателя.

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

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

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

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

Изложенные соображения и соотношения (1-5) позволяют представить в ГСК совокупность т наблюдаемых сигналов в виде вектора

>(*) =[г1(А)...2г(А)...2т(й)]Т=Ь[х;х0(А),А] + п(й), (7)

где

1-4

с

г

|Дгг(*)|

с — скорость света; (•,♦)— скалярное произведение двух векторов, указанных в круглых

скобках, определяет взаимную радиальную скорость ¿-го КА и ИИ; п(й) — /га-мерный вектор случайных дискретных последовательностей ошибок измерений частоты.

Вектор ошибок п(&) имеет гауссову ПРВ с нулевым средним значением и ковариационной

матрицей В„. Элементы главной диагонали этой матрицы есть а2{ = а2сЛ + о23, где — СКО

о

измерения несущей частоты на борту КА. Элементы вне главной диагонали равны . Взаимная ковариация М (0) - т0) • п (к)Т = 0 , поскольку при о01 « /5 и ст02 = ст03 « |Дг£ | можно пренебречь физическими предпосылками наличия статистической зависимости между координатами ИИ в априорной зоне и ошибками измерения частоты.

Алгоритм оценки координат

Задача байесовой фильтрации состояний динамической системы в постановке марковской теории нелинейной фильтрации [2, 3] при заданной функции потерь определена заданием линейных уравнений (5) и нелинейных уравнений (7). Для квадратичной функции потерь оптимальная байесова оценка х(&) текущего состояния х(&) реализуется в виде оператора апостериорного среднего по распределению вероятностей Ту(х(а) | 7,к), т.е.

оо оо

х(А)= ...]" х(Л)■ IV (х(А) | г*) <<х(Л), (8)

—оо —оо

где 2!* = |г(1), г(2— последовательность наблюдений. В силу нелинейности уравнений (7) апостериорная ПРВ в данной задаче не является гауссовой. Оператор (8) не может быть строго представлен в замкнутой форме, т.е. строго оптимальная оценка не может быть реализована на выходе дискретного линейного фильтра Калмана (ФК).

В данной задаче возможно получение квазиоптимальной текущей оценки состояния. В частности, оценка х(А) может быть реализована с помощью алгоритма обобщенного (расширенного) фильтра Калмана (РФК). Известно [2, 3], что РФК получается на основе алгоритма линейного ФК при линейной аппроксимации нелинейных функций наблюдений г (А) = Ь[х; к\ рядом Тейлора в точке известной текущей оценки состояния. Уравнения оценки и ее ковариационной матрицы имеют вид

х(*) = £"(*) +К(*)[«(*)-Г (А)]; (9)

УЖ(Л) = [Е-К(*)НЖ(А)]\^(А), (10)

где х~ (к) — экстраполированная оценка состояния; V- (А) — ковариационная матрица ошибок экстраполяции; К (к) — матричный коэффициент усиления фильтра; Нх (к) = =

= |Э/ь(«)/Ъхл __ > — матрица Якоби размерностью тпх 3. Экстраполированная оценка х (А)

[ 'ж—X (*)]

наблюдений в уравнении (9) имеет вид

Г(А) = ь[х-(А);А,ж0(А)]. (И)

Отметим два обстоятельства, связанных с реализацией алгоритма (9-11). Во первых, вычисление матрицы К (А) и оценки £~ (А) выполняются приближенно на основе линейной аппроксимации, что ведет к увеличению ошибок оценок состояния при возрастании погрешности измерения частоты. И, во-вторых, сложность выражений, определяющих функции в выражении (7), приводит к существенному увеличению объема вычислений при расчете мат-

риды Hx(fe). Желание использовать приближение 2-го порядка, как известно, требует расчета матрицы Гесса, что приводит к дополнительным громоздким расчетам.

В связи с этим целесообразно применить, развитый в работах [4, 5], метод приближенного вычисления среднего и ковариаций величин при многомерных нелинейных преобразованиях, который существенно сокращает объем вычислений. Применение этого метода к задаче фильтрации приводит к алгоритму UKF (unscented Kaiman Filter). В дословном переводе — не имеющий отношения («не пахнущий») к фильтру Калмана алгоритм.

Суть метода в том, что для вычисления (оценки) моментов нелинейной функции z = h[x] используется конечное множество точек («сигма-точки») Х„ = {хот-|, i = 0,...,2пх ; Хст е Кх, где Мх — пространство состояний и пх — его размерность. Выборка из «сигма-точек» используется для расчета всех необходимых моментов вектор-функции z = h[x] по ее выборочным значениям на множестве Х„. Расположение точек в задаче фильтрации на (А - 1)-м шаге определяется выражением

(xai(Ä-l)} = {x(Ä-l), х(А-1)±уа;(*-1)}, (12)

где ог (k -1) — i-й столбец матрицы <JVX (k -1); у = <Jnx + ^ ! ^ — составной параметр масштаба. Знак плюс в уравнении (12) принимаем для i = 1,...,пх и знак минус — для I = (пх + 1),...

..., 2пх , Квадратный корень из матрицы ^Ух (k -1) (в данном случае симметрической) следует из ее представления в виде

1=1

Для каждой «сигма-точки» (в данной задаче на этапе прогноза наблюдения) вычисляется соответствующая точка в пространстве наблюдений М2 , т.е.

z^(Ä) = h[x^(ft)]. (14)

Итоговая экстраполированная оценка наблюдений формируется по всем «сигма-точкам» (14) в виде весового среднего

z-(fe)=2X • (15)

¡=0

где сог0 = Х/(пх + X) — вес центральной (нулевой) точки хо0 (ft -1) = х (k -1); (0zi = l/[2 (nx + Л)]; i = 1 ,...,2nx — весовые коэффициенты нецентральных точек.

Для расчета матриц К(а) и V„(ä) в этом случае используется иная (эквивалентная) форма, которая связана с вычислением ковариационной матрицы Vs (й) обновляющего процесса z(&) = £z(ä) - г (ä)] и ковариационной матрицы Vxi . Указанные матрицы также вычисляются с использованием «сигма-точек». В частности:

2 пх у

vä (*) = Е С0уг • [ъ (*) - Z- (*)][z^ (*) - Г (А)] , (16)

£=0

где со у 0 = Xj[nx + X) + (l - а2 + ß2) — вес центральной (нулевой) точки; avi =cozi для всех i * 0 . Масштабный параметр X = а2 (пх +l)-nx и параметры а, ß, I обеспечивают совпадение моментов ПРВ нелинейной вектор-функции z = h[x] с моментами аппроксимирующей ее гауссовой ПРВ. Детальный анализ по этому вопросу выполнен в [4, 5].

Результаты исследования вероятностных характеристик оценок координат

Моделирование алгоритма обработки было выполнено для случая трех КА, имеющих высокоэллиптические орбиты. Наблюдаемые сигналы (7) формировались после решения уравнений (6) и расчета xQi (k) для каждого КА. Численное интегрирование уравнений (6) выполнено методом Рунге — Кутта 4-го порядка.

Координаты подспутниковых точек в момент начала интервала измерения имели следующие значения: Хх = 0°; фх = 60°N , Л2 = 120°£ ; ф2 = 60°, Х3 = 120°W; Ф3 = 60<W . Высоты КА над поверхностью Земли: fy = 4 • 104 м; йз = 4,5 ■ 104 м ; Лд = 4 • 104 м . Положение излучателя на поверхности Земли полагалось случайным в некоторой зоне S, центр которой m0 = = 84°58'W фя0 = 56°30'ЛГ]. В расчетах априорное распределение вероятностей

W"(x(0)) было выбрано равномерным. Начальные значения оценок х2 (0) = Хз0 и х3 (0) = ф50 . Дисперсии этих оценок зависят от угловых размеров зоны S, которая, в свою очередь, определяется шириной диаграмм направленности (ДН) приемных антенн, расстоянием до Земли и углами видимости на центр зоны. Угловая ширина зоны S по обеим координатам принята равной 3°.

Оценка хх (0) формируется по первым трем величинам f* = ЛГ1 [гг (l); А.а0,фз0], вычисленным на основе наблюдений z(l). Исследование вероятностных свойств величин /3* показало, что они имеют близкие к гауссовым распределения вероятностей и коррелированны между собой. Таким образом, в качестве оценки частоты fs (О) = (0) целесообразно принять максимально правдоподобную оценку f t . При выполнении численных расчетов функция правдо-

"к Г * it ♦ ^

подобия выборочного вектора fs =|/*i f$ 2 fs3 I имела гауссов вид. Максимально правдоподобная оценка fsl формировалась в виде линейной комбинацией величин f*t с коэффициентами, зависящими от элементов матрицы ковариаций. Приближенные значения ковариационных моментов вычислялись на основе Е/Т-преобразования (Unscented Transformation) [6, 7].

Исследование СКО Оф — оценок широты и а^ — СКО оценок долготы выполнено для o2fch =10; 50; 100 Гц и afs =50 Гц. Темп поступления данных в канале измерения частоты 1 с.

Результаты получены при статистическом усреднении по 200 независимым реализациям возмущений в каналах наблюдений и значений вектора х(0). На рис. 2,а,б для частоты

fs = 10 ГГц показаны зависимости СКО оценок от времени измерения.

а б

Рис. 2 — Точность оценок координат ИИ: а — СКО оценки широты; б — СКО оценки долготы

Потенциальная точность метода

Исследование потенциальной точности оценок координат предполагает вычисление дисперсий оценки (8). Отметим два способа решения этой задачи. Первый из них связан с решением уравнения Стратоновича [3, 7] для апостериорной ПРВ с последующим вычислением матрицы ковариаций оценок. Второй ориентирован на применение метода Монте-Карло для вычисления многомерного интеграла (9). В зарубежной технической литературе последний метод получения оценок (метод фильтрации) называют «particle filter» (PF — фильтр частиц) [6, 8]. Этот метод фактически следует из известного в отечественной литературе метода «важной» (существенной) выборки. Суть метода в том, что в качестве оценки апостериорной ПРВ W (х(&) | Z* j (далее используем обозначение х(&) = хк) допустимо принять эмпирическую оценку в виде

= (17)

где — нормированные весовые коэффициенты; 8(»)— многомерная дельта-функция Дира-

(0

ка; — совокупность «точечных масс» частиц в п-мерном пространстве состоянии, которые

являются независимыми выборочными значениями из существенной ПРВ | Ък). В ито-

ге в качестве оценки (8) получаем среднее арифметическое в виде

= = (18)

г=1

Известно [6, 8], что при выполнении достаточно естественных требований к существенной ПРВ (хк | г*) эмпирическая ПРВ (17) сходится почти, наверное, к истинной ПРВ. При этом (18) является асимптотически несмещенной оценкой истинного статистического среднего (8).

При условии задания существенной ПРВ в виде [6, 8]

WI(xk\Z*) = w(xk\4l1) ненормированные весовые коэффициенты связаны рекурсивным соотношением

где | х^1') — ПРВ, определяющая правдоподобие наблюдений для частицы .

В таблице приведены результаты фильтрации координат для = 10 Гц и о/5 = 50 Гц при

использовании алгоритмов ЩУ? и РЕ. Показаны СКО оценок координат ИИ при усреднении по 100 независимым реализациям на тридцатом временном шаге (Л = 30 ).

Алгоритм о^, угл. мин а,- , угл. мин Ао Примечание

UKF 1,739 2,644 а = 0,1; р = 2; / = 1

PF 2,631 3,383 Число «частиц» N = 24 000

PF 1,827 2,664 Число «частиц» N = 98 000

PF 1,782 2,641 Число «частиц» N = 250 000

Результаты моделирования показывают сходимость расчетных значений СКО оценок координат при увеличении количества частиц в алгоритме PF. Отметим, что решение задачи с двумя КА дает в конечной точке интервала наблюдения (А = 100). СКО оценок для обеих координат примерно 24' и одного КА — около 50'.

Заключение

Предложенный в работе алгоритм обработки измерений текущей частоты радиосигналов, принятых на борту одного или нескольких КА, находящихся на высокоэллиптических орбитах, обеспечивает формирование устойчивых оценок координат наземного ИРИ с неизвестной частотой излучения. Точность оценок координат при СКО измерения частоты ofch = 10 Гц и Ofs =50 Гц фактически соответствует потенциально достижимой точности байесовых оценок при квадратичной функции потерь.

Выполнение всех вычислений на текущем временном шаге требует 1,5 мс при использовании ЭВМ с процессором Athlon 64 3000+.

Литература

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

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

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

4. Julier, S.J. The Scaled Unscented Transformation / S.J. Julier // Proceedings of the American Control Conference. - V. 6, Anchorage, AK, USA, May 200. - P. 4555-4559.

5. Julier, S.J. A new approach for filtering nonlinear systems / S.J. Julier, Uhlmann and Durrant-Whyte // Proceedings of the American Control Conference. - 1995. - P. 1628-1632.

6. Geweke, J. Bayesian Inference in Econometrics Models using Monte Carlo Integration / J. Geweke. - Econometrica. - 1989. - V. 57.

7. Стратонович P.JI. Условные марковские процессы и их применение к теории оптимального управления / P.JI. Стратонович. - М. : МГУ, 1966.

8. Doucet, A. Monte Carlo Methods for Bayesian Estimation of Hidden Markov Models / A. Doucet // Application to Radiation Signals, PhD. Thesis, Univ. Paris-Sud, Orsay, 1997.

9. Сильвестров С.Д. Точность измерения параметров движения космических аппаратов / С.Д. Сильвестров [и др.] ; под ред. С.Д. Сильвестрова. - М. : Сов. Радио, 1970.

Савин Александр Александрович Аспирант каф. радиотехнических систем ТУСУРа Телефон: (3822) 41 38 92 Эл. почта: [email protected]

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

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

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

Эл. почта: [email protected]

Savin А.А., Tislenko V.I.

A quasi-optimal estimation of the coordinates of a ground radiation source in space-based system measuring radio signals frequency

A relative openness of space links used for information transferring gives a real possibility

for various malefactors to have an unauthorized access to these links. A radical way of struggling against

such «pirates» assumes determination of a ground radiation source coordinates with the use of signals

which can be received by one or more satellites. Calculation results of root-mean-square errors (RMS)

of joint estimates for a ground radiation source's coordinates (latitude, longitude) and frequency

are presented in the paper. The estimates are generated at the output of a quasi-optimal filter,

which input signals are measurements of HF carriers of the ground radiation source received

by onboard receivers of one or more satellites at high-elliptic orbits. RMS values of the quasi-optimal

estimates are compared with the values of optimal Bayes estimates provided quadratic loss.