Устойчивый алгоритм предварительной обработки измерений псевдодальностей в системе ГЛОНАСС Текст научной статьи по специальности «Математика»

Научный вестник НГТУ. - 2009. - № 3(36)

ОБРАБОТКА ИНФОРМАЦИИ

УДК 629.78: 517.518.8

Устойчивый алгоритм предварительной обработки измерений псевдодальностей в системе ГЛОНАСС*

Ю.Е. ВОСКОБОЙНИКОВ, А.С. ТОЛСТИКОВ

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

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

ВВЕДЕНИЕ

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

Происхождение выбросов или пачек выбросов обычно связывают с отражениями навигационных сигналов от расположенных вблизи приемной антенны предметов и объектов, что создает эффект удлинения радиотрассы. Пропуски данных возникают вследствие потери синхронизации приемной аппаратуры на отдельных интервалах времени. Зачастую дефекты траекторных измерений (т. е. аномальные измерения и пропуски данных), особенно при работе приемной аппаратуры в городских условиях, составляют 15-25 % от общего объема данных. Некоторые задачи (например, в практике геодезических определений), допускают отбраковку (удаление из дальнейшей обработки) всех данных, содержащих дефектные измерения. Однако существуют задачи координатно-временных определений, в которых такая отбраковка непозволительна. Так, в задачах синхронизации группы пространственно -разнесенных часов отбраковка результатов измерений в одном из пунктов разрушает весь сеанс синхронизации. При решении задач формирования эфемеридно-временного обеспечения (ЭВО) спутниковых навигационных систем потеря данных траекторных измерений на отдельных дугах орбиты внутри мерного интервала существенно снижает качество ЭВО.

Заметим, что измерение геометрических дальностей от навигационного спутника до приемной аппаратуры беззапросными методами сводится к измерениям интервала времени, необходимого для прохождения навигационного сигнала от передающей антенны спутника до антенны приемной аппаратуры. Этот интервал времени, выраженный в единицах длины, и называется псевдодальностью. Поскольку аномальные измерения имеют небольшой уровень (порядка 100 м) по сравнению с измеряемыми дальностями (порядка 19000 — 25 000 км ), то переходят к рассмотрению движения навигационного спутника относительно некоторой опор-

* Статья получена 10 февраля 2009 г.

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

■А /

1, "Уй тп 11 т

1

0 1 2 3 4 Часы

Рис. 1. Исходные данные и результат фильтрации медианным фильтром

Таким образом, описанная модель измерений псевдодальностей требует решения двух задач предварительной обработки:

• фильтрации(удаления) аномальных измерений,

• гладкого восполнения пропущенных измерений.

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

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

1. НЕЛИНЕЙНАЯ ЛОКАЛЬНАЯ ФИЛЬТРАЦИЯ ИЗМЕРЕНИЙ ПСЕВДОДАЛЬНОСТЕЙ

Аномальные измерения псевдодальностей можно рассматривать как импульсивные шумы, амплитуда которых в несколько раз может превосходить значения соседних измерений. Для фильтрации таких шумов целесообразно использовать нелинейные локальные фильтры, в частности медианный фильтр (МФ) [1, 2]. Для пояснения принципа работы этого фильтра предположим, что даны значения некоторого зашумленного одномерного сигнала:

/у = /у + 1 < у < И, где у - случайный шум. Тогда выходной сигнал /уМФ медианного

фильтра определяется соотношением

/ ,МФ = ше^(, 1_ь+1,.., / ,,...,}, +ь), (1)

где medL () - функция, вычисляющая медиану из значений, указанных в скобках. Таким образом, fjМФ равна медиане выборки из 2L+1 значений fj , сформированной из L значений,

лежащих левее точки j, L значений правее точки j и значения fj . Величину 2L+1 можно

интерпретировать как размер апертуры медианного фильтра. «Перемещая» апертуру фильтра

f f МФ

по всем значениям f , вычисляют значения f для всего исходного сигнала.

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

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

• средние значения и дисперсия остаточного шума определяется соответствующими числовыми характеристиками шума «нормальных» измерений.

Известно, что медианный фильтр хорошо фильтрует импульсный шум, но плохо - низкоамплитудный. В нашем случае это шумы нормальных измерений. Для фильтрации таких шумов в [3] предложен класс комбинированных фильтров, представляющих собой последовательное соединение медианного фильтра (удаление импульсного шума) и некоторой модификации фильтра скользящего среднего (удаление низкоамплитудного шума с сохранением контрастных деталей, присутствующих в «точном» сигнале fj ). Под контрастными деталями будем понимать резкое изменение амплитуды сигнала, что обусловливает большие производные сигнала. Работу такого комбинированного фильтра (КФ) «медианный фильтр + интервальный фильтр скользящего среднего» можно представить следующими шагами, выполняемыми для j е [1,..., N]:

Шаг 1. Строится оценка:

fjМФ = med L (Г}_ьfjfJ+L). (2)

Шаг 2. Строится оценка:

fjКФ = ave% (f МФ : j - K < i < j + K, |f МФ - fjМФ | < р), (3)

где averK () - операция усреднения, L, K - размеры апертур фильтров, причем K > L . Заме-

~ /-МФ

тим, что на шаге 2 усредняются только те значения f , которые попали в интервал

Ар (j) = |/уМФ-Р, fjM° + р] (4)

Такое интервальное усреднение предотвращает сглаживание контрастных деталей сигнала fj и поэтому выбор величины Р играет существенную роль в работе комбинированного фильтра. Рассмотрим задание Р по дисперсии CTj2 . Предположим, что шум j нормальных измерений имеет дисперсию ст2 . Тогда (используя известное правило «трех сигм»)

Р = 3а§. (5)

2

Соотношение (5) требует задания дисперсии ст^ . В большинстве случаев обработки реальных измерений эта величина неизвестна и для ее оценивания рекомендуется использовать приводимую ниже методику.

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

2. ВОСПОЛНЕНИЕ ПРОПУЩЕННЫХ ИЗМЕРЕНИЙ ПСЕВДОДАЛЬНОСТЕЙ

Восполнение пропущенных измерений можно осуществить, используя интерполяционный многочлен, построенный по соседним (к пропущенным данным) измерениям. В качестве такого интерполянта будем использовать кубический сглаживающий сплайн (СКС) по следующим причинам [4,5].

1. Кубический сплайн дает интерполянт с непрерывной второй производной на всем интервале значений псевдодальностей, что позволит не только восполнить пропущенные измерения, но при необходимости (при дальнейшей обработке) вычислить производные до второго порядка включительно.

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

Пусть функция ф(т) задана на интервале [т^тт ] таблицей {т;,фг} , г = 1,2,...,т , где

фг = ф(т;) + "Пг, ■ = 1,2,...,т , (6)

2

- случайная погрешность с нулевым средним и дисперсией ст^. . Сглаживающий кубический сплайн ¿а (т) на каждом интервале [т;, тг+1 ] является кубическим многочленом вида

¿а (т) = а. + Ьг (Т _ Т; ) + Сг (Т _ Т; )2 + ¿г (Т _ Т ) . (7)

Для однозначного вычисления коэффициентов сплайна а1, , с(, дополнительно задаются краевые условия. Для рассматриваемой задачи восполнения пропущенных измерений примем так называемые «естественные» краевые условия вида

¿а(т )=¿ж )=о. (8)

Можно показать [4,5], что СКС с условиями (8) доставляет минимум сглаживающему функционалу

Тп 2 т 2

Ма [5] = а | (Лх))2 ¿х + X р_ (Фг _ ))2 , (9)

Т1 '■=!

где а > 0 - параметр сглаживания, р1 > 0 - весовые множители, пропорциональные точности измерения ф; (в частности, р1 □ ст2. , г = 1, 2,..., т).

При заданном параметре сглаживания а коэффициенты а;, Ь;, с;, определяются через вектор решения системы линейных алгебраических уравнений с пятидиагональной матрицей. Элементы этой матрицы зависят от а . При а ^ 0 сглаживающий сплайн становится интерполяционным и, следовательно, его значения содержат погрешности "П.. При а ^ ж сплайн переходит в прямую линию (вторая производная равна 0 на всем интервале [т^ тп ]).

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

задачи восполнения в качестве такого функционала примем среднеквадратическую ошибку (СКО) сглаживания:

т 2

Д($х) = М К^Т, )"ф(хг ))

_г=1 _

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

Оценивание аор( на основе критерия оптимальности сглаживающего сплайна. Алгоритм подробно описан в [5, 6]. Поэтому здесь только заметим, что в качестве оценки для аор1 принимается значение (обозначим а^ ), при котором величина

р„, (а) = £

,=1 аЛ1

попадает в интервал [Ут,0025, Ут 0975 ] , где Ут,0025 , Ут 0975 - квантили х2 -распределения

с т степенями свободы уровней 0.025, 0.975 соответственно. Доказывается существование

2

а^ при любых дисперсиях ст^., сходимость СКС, построенного при а = а^ к интерполяционному при ст2. ^ 0 [5]. В вычислительных экспериментах с различными функциями ф(т) показано, что увеличение СКО по сравнению с аор1 не превышает 10 —15 % .

Оценивание аор( на основе метода Ь -кривой. Метод Ь -кривой используется для выбора параметра регуляризации при решении некорректно-поставленных задач (например, см. [7]). Суть этого метода для оценивания аор1 сглаживающего сплайна заключается в следующем. Введем величины

т 1 2 Тт 2

р(а) = Х Р— (ф, — 5г (т,)) , Т(а) =К^ (Т)) .

Параметрическая функция с координатами (р(а), у(а)) называется Ь -кривой. В качестве оценки аор1 принимается значение а (обозначим аь ), при котором Ь -кривая имеет максимальную кривизну.

Численная реализация этого метода, требующая вычисления производных р'(а) , у'(а), приведена в [7]. Здесь же выполнено исследование точности оценки аь. Следует отметить, что вычисление аЬ требует на порядок больше вычислительных операций и СКО сплайна, построенного при а = аь , большей, чем при а^ , но зато в методе Ь -кривой не требуется

2

задания дисперсий ст^ .

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

Шаг 1. Формируется два массива: массив |фу | из отфильтрованных непропущенных измерений /МФ (или /¡КФ) и массив |ту | из аргументов значений фу . Длина т массивов

ху , фу меньше (при наличии пропусков) длины п и на практике доля пропущенных измерений (п — т)/п может достигать 15 — 20 % .

Шаг 2. По таблице |ху, фу |, у = 1, 2,..., т строится СКС, параметр сглаживания которого равен вычисленному а^ (или а ь ) и определяются коэффициенты а у , Ьу , Су , ёу .

Шаг 3. Определяется значение хр из массива |ху |, правее которого находятся аргументы пропущенных измерений.

Шаг 4. Вычисляются значения сплайна £"а (х) для значений пропущенных измерений (см. (7))

8а(1, ) = ар + Ьр —хр ) + Ср —хр ) + ёр —хр ) . <10)

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

Оценивание дисперсии. Для выбора параметра сглаживания на основе критерия опти-

2

мальности требуется задание дисперсий ст^. . Ранее было показано, что остаточный шум ме-

2 ~ дианной фильтрации имеет дисперсию ст^ погрешностей исходных данных . Поэтому при

2 2 2 выборе параметра а полагаем, что ст^ = ст^ . Если дисперсия ст^ неизвестна, то для оценивания дисперсии предлагается алгоритм [7, стр. 102-103], который можно представить следующими шагами.

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

Шаг 2. Из проекций вектора 2 формируется периодическая (с периодом Nр) последовательность по правилу:

.(у ) =

= у+1, у = 0,..., N — 1, р (у) 1 о, у = N,..., Nv — 1.

р

Величина периода Nр > N обусловлена дальнейшим применением алгоритма быстрого преобразования Фурье и очень часто является степенью 2, т. е. 2р . Тогда в качестве Nр принимается наименьшее 2р > N .

Шаг 3. Вычисляются коэффициенты Xр (I) дискретного преобразования Фурье (ДПФ)

последовательности (у с использованием алгоритма быстрого преобразования Фурье

1 NP —1 1 2то 1

2р (1 ) = N I 2р (у)ехр 1— • у|, I = 0,...,^ — 1.

Шаг 4. Вычисляется оценка

2 N1 ~ / ч|2

^) = М^^р (Мр'2+')1 .

Эта оценка основана на следующем свойстве ДПФ. Величины р (I)| симметричны от. |2

носительно точки Nр /2 и в некоторой окрестности этой точки значения \2р уIV обусловлены только случайными погрешностями. Величина задается равной (0.1...0.15)Nр . При-

2

веденный алгоритм исследован в [7, стр. 103] и позволяет оценивать дисперсию ст^ с относительной погрешностью 5 —10 % . Этого достаточно для надежного определения аор1 первым алгоритмом.

3. РЕЗУЛЬТАТЫ ОБРАБОТКИ ЭКСПЕРИМЕНТАЛЬНЫХ ДАННЫХ

Для иллюстрации работоспособности предложенного алгоритма первичной обработки приведем некоторые результаты следующего эксперимента. В качестве исходных данных были взяты реальные измерения псевдодальностей. Количество измерений п = 723. На рис. 1 эти измерения показаны штриховой линией. Видно наличие аномальных измерений и одного интервала пропущенных измерений, количество которых равно 131. Приведенные измерения были подвергнуты медианной фильтрации (1) с величиной Ь = 15 . Отфильтрованные значения приведены на рис. 1 (сплошная линия). Видно, что все аномальные измерения, присутствующие в исходных данных, успешно удалены. При этом уменьшился и уровень шума нормальных измерений псевдодальностей.

Затем по отфильтрованным значениям псевдодальностей был построен сглаживающий ку-

2

бический сплайн с параметром сглаживания а = а^ . Так как дисперсия ст^ была неизвестна,

то по первым 320 измерениям была вычислена оценка (Np = 512, Ь = 50 ), которая использовалась при вычислении а-^ . После этого по формуле (10) были определены значения сплайна для аргументов пропущенных измерений. На рис. 2 сплошной кривой показаны значения восполненных значений псевдодальностей, а точками - отфильтрованные значения псевдодальностей, по которым строился сглаживающий кубический сплайн. Видно, что применение СКС позволило не только восполнить пропущенные измерения, но и уменьшить уровень погрешностей непропущенных измерений.

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

200

100

-100

-200

/

0 12 3

Часы

Рис. 2. Результаты восполнения пропущенных измерений

эксперимент. С помощью программного имитатора измерительной информации [8] были получены фрагменты траектории навигационного спутника ГЛОНАСС, по которым вычислялись псевдодальности. Затем эти значения исказили случайными погрешностями, включая импульсные (моделирование аномальных измерений), и 18 % измерений удалили. К полученным таким образом данным был применен предлагаемый алгоритм обработки. Сравнение исходных «точных» значений псевдодальностей с восстановленными значениями показало высокую точность восполнения - погрешность восполнения не превысила 20 см.

ЗАКЛЮЧЕНИЕ

Предложенный в работе подход позволил успешно решить сложные задачи первичной обработки результатов измерений псевдодальностей. При этом потребовался минимум априорной информации об обрабатываемом сигнале (предполагалось, что сигнал на каждом из интервалов дискретизации [i7-, t/+1 ] хорошо аппроксимируется кубическим полиномом). Алгоритм может

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

СПИСОК ЛИТЕРАТУРЫ

[1] Глобальная навигационная спутниковая система ГЛОНАСС: Интерфейсный контрольный документ (редакция 5.0). - М.: 2002.

[2] Претт У. Цифровая обработка изображений. - М.: Мир, 1982. - Кн. 2.

[3] Воскобойников Ю.Е. Нелинейные комбинированные алгоритмы фильтрации зашумленных сигналов и изображения // Автометрия. - 1990. - № 1. - С. 21-28.

[4] Завьялов Ю.С., Б.И.Квасов, В.ЛМирошниченко. Методы сплайн-функций. - Новосибирск: Наука, 1980.

[5] Воскобойников Ю.Е., Преображенский Н.Г., Седельников А.И. Математическая обработка эксперимента в молекулярной газодинамике. - Новосибирск: Наука, 1984. -

[6] Воскобойников Ю.Е. Оценивание оптимального параметра регуляризирующих алгоритмов восстановления изображения // Автометрия. - 1995. - № 3. - С. 68-79.

[7] Воскобойников Ю.Е. Устойчивые методы и алгоритмы параметрической идентификации: монография. -Новосибирск: Изд-во НГАСУ. 2006 (электронная версия http://old.sibstrin.ru/prikl/monogr07.html)

[8] Владимиров В.М., Гречкосеев А.К., Толстиков А.С. Имитатор измерительной информации для отработки эфемеридно-временного обеспечения космической навигационной системы ГЛОНАСС // Измерительная техника. -2004. - № 8. - С. 12-14.

Воскобойников Юрий Евгеньевич, доктор физико-математических наук, профессор, заслуженный работник высшей школы РФ, заведующий кафедрой прикладной математики Новосибирского государственного архитектурно-строительного университета, профессор кафедры автоматики Новосибирского государственного технического университета. Академик Международной академии информатизации, Российской академии естествознания, член-корреспондент СО АН ВШ. Основное направление научных исследований - методы и алгоритмы решения некорректных задач интерпретации экспериментальных данных, нелинейные алгоритмы фильтрации, вычислительная томография. Имеет более 240 публикаций, в числе которых три монографии и 14 учебных пособий.

Толстиков Александр Сергеевич, кандидат технических наук, начальник отдела «Измерения времени, частоты и параметров вращения Земли» Сибирского НИИ метрологии, доцент кафедры метрологии, сертификации и стандартизации Сибирской государственной геодезической академии. Член-корреспондент Российской метрологической академии. Основные направления научных исследований - метрологическое обеспечение спутниковых навигационных технологий, алгоритмы обработки беззапросных траекторных измерений .Имеет более 80 публикаций.

Yu. E. Voskoboinikov, A.S.Tolstikov

Steady algorithm preliminary processing of pseudorange measurements in GLONASS

The offered algorithm of processing of pseudorange measurements solves two problems of primary processing: a filtration of abnormal measurements and completion of the missed pseudorange measurements. The first problem is solved with use of the nonlinear local filters, the second - with application of smoothing cubic splines. The results of processing of real measurements and the data received on the program simulator of measuring information GLONASS, have shown efficiency of the offered algorithm.

Key words: trajektory measurements, pseudorange measurements, filtration of abnormal measurements, completion of the missed measurements, the nonlinear filtration, a smoothing cubic spline.