Вычисление положения искусственного спутника Земли при ограничениях на ошибки измерений Текст научной статьи по специальности «Математика»

E. A. Popov, M. E. Semenkina, L. V. Lipinskiy

DECISION MAKING WITH INTELLIGENT INFORMATION TECHNOLOGY ENSEMBLE

Genetic programming algorithm for neural network automatic design is suggested. Ensemble member competence estimations, based on decision making with intelligent information technologies, is proposed. The approach effectiveness is approved on qualificatory and real-world problems.

Keywords: intelligent information technologies, genetic programming algorithm, group decision making.

© Попов Е. А., Семенкина М. Е., Липинский Л. В., 2012

УДК 517.977.1

А. Н. Рогалев, А. А. Рогалев

ВЫЧИСЛЕНИЕ ПОЛОЖЕНИЯ ИСКУССТВЕННОГО СПУТНИКА ЗЕМЛИ ПРИ ОГРАНИЧЕНИЯХ НА ОШИБКИ ИЗМЕРЕНИЙ

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

Ключевые слова: метод наименьших квадратов, минимаксный метод, фильтрация данных.

В задачах оценки параметров объектов астрометрии и космической геодезии необходимо анализировать экспериментальные данные, полученные с помощью измерений [1-4]. Поскольку позиционные наблюдения, как правило, являются косвенными, то результаты (данные) часто представляются алгебраическими моделями с неизвестными коэффициентами, которые являются искомыми параметрами. При оценивании этих параметров в условиях, когда данные наблюдений получены с ошибками, применяются подходы, основанные на методе наименьших квадратов [4-6] и минимаксном методе [4; 7-10].

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

пределение, которое хорошо аппроксимируется равномерным законом. Эта постановка не соответствует предположениям, при которых доказаны оптимальные свойства оценок МНК.

Проводится сравнительный анализ нескольких вычислительных схем, построенных на основе применения МНК и минимаксного метода (МНМ) в задаче вычисления положения искусственного спутника Земли при ограничениях на ошибки измерений. Следует отметить, что значимость этих результатов определяется тем, что не всегда ясно, когда эффективен тот или иной подход к решению задачи: классический, при котором характеристики случайных ошибок исходных данных основываются на методе максимального правдоподобия, и неклассический (минимаксный), характеризуемый заданием лишь некоторых множеств, которым принадлежат ошибки или их характеристики. Во втором случае интересно также сравнение возможностей минимаксного и гарантированного подходов [11; 12] к оценке решений управляемых систем.

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

Предположим, что движение центра масс КА в околоземном пространстве с достаточной степенью точности и адекватности описывается следующей

системой уравнении в относительной геоцентрической системе координат OAXYZ :

dx _

1ґ = ^ ’ dy dt у ’

dv,

л = 1°з + Г^1x+’

dz

dt

■ = V,

dt

dt

к2 =7^

(1)

где

Я 3 (я

к1 = аоо з а зо I

к =- Я _ 3 (я

кз = аоо г з азо I

Л

5 % _ 1

Г

V ' У

V .з

Л

5 % _ -

V Г /

V ' У

аоо = '

-азо =

з Г = 4 хз + уз + ,з.

„„ Я Я3

Введем обозначения: п0,п2 - коэффициенты разложения потенциала ускорения силы притяжения Земли в ряд по сферическим функциям Лежандра; Я - средний радиус Земли; - угловая скорость

вращения Земли.

Пусть имеются также расчетные значения

Хз > Уо> ,v

(з)

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

Расчетные сведения (з) могут быть заданы из априорных сведений или найдены как оценки в результате предварительного определения орбиты КА по достаточному количеству измерений [з; 4; 7].

В качестве измеряемых параметров рассмотрим наклонные дальности, измеряемые с трех радиолокационных станций, расположенных на поверхности Земли в точках Р1 (х1,у1,,х), Рз (хз,уз,,2),Р3(х3,у3,,3) с известными координатами. Уравнения для определяемых параметров имеют вид

Р1 = 4(х- х1)2 + (-У->1)2 + (2 -2г)2;

Р2 = 4(Х - Х2)2 + (У - У2)2 + (2 - г2)2 ;

Рз = 4(х - хз)2 +(У - Уз)2 +(2 - 23)2.

Предположим, что измерения всеми станциями производятся синхронно в дискретные моменты времени ,1 = 1,2,..., N, с постоянным шагом № по времени. Результаты измерений определяются по формулам

Рз

Рз

(^) = 4 (х _ ^)з + (у _ ^)з + (, _ ^)з +Дрі,г;

) = 4(х _ хз)з + (у _ Уз)з + ^ _ ,з)з + ДР з,г;

(^ ) =у] (х _ Х3 )з + (у _ У3 )з + (, _ ,3 )з + Дрз,і,

Движение спутника будет определено в пространстве, если известны плоскость, в которой лежит его орбита, размеры и форма этой орбиты, ее ориентация в пространстве и момент времени, в который спутник находится в определенной точке орбиты. Чтобы охарактеризовать ориентацию орбиты в пространстве, нужно прежде всего задать базовую систему координат, начало которой совпадает с фокусом орбиты О (т. е. с центром Земли). За основную плоскость ОХУ, относительно которой определяется положение орбиты спутника, принимается плоскость экватора. Ось ОХ пересекает экватор в точке с долготой 0°, соответствующей Гринвичскому меридиану.

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

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

Современные средства обработки косвенных измерений позволяют получать результаты с очень высокой точностью до 10-6. В работе [2] отмечено следующее: чтобы не потерять эту точность, нужно иметь возможность вычислять пространственные положения ИСЗ с методической точностью, по крайней мере в три раза превышающей точность наблюдений.

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

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

где Др1 і, Дрз і, Др3 і _ ошибки измерений.

4 =^,.(01,0з,..., 0„) + Ді, і = 1,з,

(3)

іоо

где ё, - результат /-го измерения, содержащий случайную ошибку; ё, - вектор размерности N, чаще всего делается несколько косвенных измерений; Ь (01,..., 0 N ) - некоторая известная векторная функция от неизвестных параметров; разность

Д, = 4 -^ ^са^ 02,са1с , ..., ^сак^ - вектор невязки

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

Оценками, получаемыми методом наименьших квадратов (МНК-оценками), называется совокупность значений 01, 02 , ..., 0N , обращающих в минимум взвешенную сумму квадратов невязок

Ф(01, 02, ..., 0N ) = ^Срг [ -[, (01, 02 , ..., 0N )] .

,=1

Для прямых измерений

ф(01, 02,..., 0N ) = ]Гр, [ -0, ]2.

, =1

Если функции Ь, (01 , 02 ,..., 0N ) нелинейны, то производится их разложение в ряд Тейлора и отбрасывание членов второго и более высоких порядков. Тогда уравнения, входящие в фундаментальную систему, можно записать в виде

d1 =£,

..0У)+1 %

Д0 + Д

0=0(о)

Н, К - ортогональные

матрицы,

Я - блочно-

диагональная матрица.

Применение метода МНМ основано на минимизации суммы модулей невязок вектора измерений.

Пусть требуется найти вектор ©, который бы обращал в минимум выражение

Н (©) = Н (0!, 02,..., 0N ) =

П . (4)

= Ё _^ (01, 0з , ..., 0 N ^

,=1

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

Если провести линеаризацию фундаментальной системы уравнений, то вместо (4) следует минимизировать выражение

Ф(,) = Х Рг

і =1

Щ _Ё V.

} =1

(5)

гдЄ ^ = 0 _0

(о)

Ыг = йг -£

(о),

Ь„ =

50,

0=0(о)

Система уравнений для определения значений 01,0з,...,0N, содержащая N уравнений относительно N неизвестных, называется системой нормальных уравнений. Возможны различные варианты решения нормальной системы уравнений при ее линеаризации. Был выбран подход, хорошо изложенный в [13]. Основная линейная задача наименьших квадратов формулируется следующим образом.

Пусть даны действительная (т х п)-матрица А ранга к < тіп(т, п) и действительный т-вектор Ь . Найти действительный п-вектор хо, минимизирующий евклидову длину вектора Ах - Ь.

Решение основано на представлении (т х п)-мат-

рицы А в виде произведения А = НЯКТ, где

Для минимизации выражения (5) можно применить два подхода. Первый из них связан с использованием идей линейного программирования, а второй -с использованием вариационно-взвешенных квадратических приближений [10], являющихся обобщением вычислительной схемы метода наименьших квадратов. Алгоритмы минимизации суммы модулей невязок, построенные на базе идей линейного и кусочнолинейного программирования, имеют одну привлекательную черту: они позволяют легко наложить на искомые оценки некоторые ограничения в виде неравенств, если заранее известно, что оценки не должны выходить за определенные пределы. Задача минимизации суммы модулей невязок представляет собой одну из самых простых задач кусочно-линейного программирования, когда отсутствуют какие-либо ограничения типа gi ^, г2,..., ) < ё, на переменные 2,.

Эта задача эквивалентна следующей задаче линейного программирования:

Фі =Ё Рг (иг + V ), і=1

+ иі _Vi = ° і = 1,3,

(6)

. = 1,з,..., N; и > о, V > о,

,Ф - ,(2) / =

] ] ]

і = 1,з,..., п; ,(1) > о, ,(2) > о,. = 1,з,..., N.

(7)

_ ’ ./

Доказательство эквивалентности следует из того, что щ, V, одновременно не могут быть положительны, так как иначе выражение (6) могло бы быть уменьшено на величину 5, = 2р, шт(и,, V,) без какого бы то ни было нарушения ограничений, входящих в (7). Поэтому

мі _Е .

.=1

= \иі _ vi\ = и + vi

Ю1

и решение системы неравенств (12), минимизирую-

п

щее Ф, = ^ р1 (и, + vi), одновременно минимизирует

,=1

интересующую нас сумму модулей невязок.

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

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

ем л (8)

,=1 р, \ёг -Ь (V)|

Очевидно, что

Н (©) = ]Гр, |4 -Ь (01, 02,..., 0 N )| = е(©, ©) .

, =1

Пусть ©(0) - некоторое начальное приближение. Подставим в (8) вместо вектора V фиксированный вектор ©(0). Тогда 0(©,©(0)) будет представлять собой квадратическую форму относительно вектора ©. Вектор ©(1), обращающий в минимум форму е(©, ©(0)), можно найти, используя вычислительную схему метода наименьших квадратов. Далее рассмотрим квадратическую форму е(©, ©(1)) и найдем вектор ©(2), обращающий в минимум эту форму и т. д. Можно показать, что полученная таким образом последовательность векторов ©(0),©(1),..., ©('?'1,... схо-*

дится к вектору © , минимизирующему сумму модулей. Подстановка в выражение (8) вместо вектора V различных значений ©(0), ©(1),..., ©(х),... приводит фактически к тому, что на каждом этапе вычислений векторов производится корректировка (вариация) весов, так что ■?) =|^р\й1 -Ь, (©(*))|^ . Поэтому рассматриваемый метод получил название метода вариационно-взвешенных квадратических приближений.

Известно, что если функции Ь (©) нелинейны, то *

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

©(0), ©(1),..., ©(х),... потребовалось бы несколько ите-

раций, т. е. в общей сложности время счета должно было возрасти в несколько раз.

Но оказывается, что вовсе необязательно минимизировать форму 0(©, ©(г^) по © . Достаточно получить хотя бы один вектор ©(г+:) такой, что е(©(г-1, ©(г-1) > е(©(г+1-1, ©(г-1), и процедуру отыскания последовательности векторов ©(г+2), ©(г+3),..., сходящихся к искомому вектору ©п, можно продолжить.

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

Для сравнения алгоритмов этих подходов моделируются различные варианты получения радиолокационных измерений в задаче вычисления положения искусственного спутника Земли при ограничениях на ошибки измерений, отличающихся между собой числом проходов космического объекта через зону одного или нескольких радаров (см. рисунок). Модель движения задана системой (1) [з; 14].

Вычисление положения ИЗС

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

Большая полнота сведений о погрешностях исходных данных, используемых при построении алгоритма

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

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

Таблица 1

Улучшение оценок параметров при сложных распределениях ошибок за счет использования МНМ

Параметры Относительное улучшение точности оценки МНМ по сравнению с МНК для сложного распределения, %

X 6,23

Y 78,15

Z 81,46

Vx 40,96

Vy 76,55

Vz 96,34

Таблица 2

Точность оценок параметров движения при равномерном распределении ошибок измерений в пределах 50 м

Параметры Итоговая точность оценки, м

Метод

МНК МНМ

X -0,433 -1,42

Y 1,343 3,72

Z -2,917 -6,1

Vx -0,009 -0,008

V -0,058 -0,084

Vz 0,05 0,077

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

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

Библиографические ссылки

1. Бордовицына Т. В. Современные численные методы в задачах небесной механики. М. : Наука, 1984.

2. Бордовицына Т. В., Авдюшев В. А. Теория движения искусственных спутников Земли. Аналитические и численные методы : учеб. пособие. Томск : Изд-во Том. ун-та, 2007.

3. Дубошин Г. Н. Небесная механика. Основные задачи и методы. М. : Наука, 1968.

4. Эльяберг П. Е. Определение движения по результатам измерений. М. : Наука, 1976.

5. Линник Ю. В. Метод наименьших квадратов и основы математико-статистической теории обработки наблюдений. М. : Физматгиз, 1962.

6. Губанов В. С. Обобщенный метод наименьших квадратов. Теория и применение в астрометрии. СПб. : Наука, 1997.

7. Самотохин А. С., Хуторовский З. Н. Определение прогнозного положения ИСЗ при ограниченных ошибках измерений // Космич. исслед. 2011. Т. 49. № 6. С. 526-537.

8. Черницов А. М., Тамаров В. А. Определение допустимого уровня систематических ошибок наблюдения при построении областей возможных движений малых тел // Изв. вузов. Физика. Прил. Небесная механика и прикл. астрономия. Томск, 2006. Т. 49, № 2. С. 75-81.

9. Лидов М. Л. Минимаксные методы оценивания

[Электронный ресурс] // Препр. ИПМ им. М. В. Келдыша. 2010. № 71. URL:

http://library.keldysh.ru/preprint.asp (дата обращения: 12.12.2012).

10. Мудров В. И., Кушко В. Л. Методы обработки измерений. Квазиправдоподобные оценки. М. : Радио и связь, 1983.

11. Rogalyov A. N. Computation of reachable sets guaranteed bounds // Proc. of the IASTED Intern. Conf. on Automation, Control, and Information Technology -Control, Diagnostics, and Automation. Calgary, 2010. P. 132-139.

12. Рогалев А. Н. Вопросы реализации гарантированных методов включения выживающих траекторий управляемых систем // Вестник СибГАУ. 2011. Вып. 2(35). С. 54-58.

13. Лоусон Ч., Хенсон Р. Численное решение задач наименьших квадратов. М. : Наука, 1986.

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

A. N. Rogalyov, A. A. Rogalyov

NUMERICAL COMPUTATIONS FOR ARTIFICIAL EARTH SATELLITE POSITION UNDER RESTRICTIONS ON MEASUREMENTS ERRORS

The application of the least squares and minimax methods for near-earth space object position identification under known restrictions on the measurement errors is considered in this article.

The computational schemes being used for a space object trajectories parameters evaluating under different measurement errors are studied.

The given algorithms accuracy characteristics under different measurement errors distributions are determined. The results of the computations are presented.

Keywords: method of the least squares, minimax method, motion trajectory identification, data filtering.

© PoraneB A. H., PoraneB A. A., 2012

УДК 517.95

С. И. Сенашов

О СИСТЕМАХ ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ С ДВУМЯ ХАРАКТЕРИСТИКАМИ

Показано, что система т (т > 2) квазилинейных дифференциальных уравнений от двух переменных, имеющая две характеристики, сводится к системе двух квазилинейных уравнений и системе т - 2 линейных уравнения.

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

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

а.у (и)дхи] + Ьу (и)дуи3 = 0, і, ] = 1,..., т, (1)

где и = (и1,..., ит).

Пусть система (1) имеет только две характеристики

± = лй, іу = БЙ,

dx іх

а также т инвариантов Римана I1,..., Ет.

В этом случае система (1) может быть записана в виде

+ Л(^,..., 1т) = 0,

дх ду

+ Б(^,..., Ет) = 0,

дх ду

+ С(^,..., 1т) — = 0, і = 3,..., т,

дх ду

где С равно либо А, либо В.

Для простоты рассмотрим случай т = 4, который часто встречается, например, в теории пластичности. Имеем

+ Л(^,..., |4) = 0,

дх ду

^+E(I1,

дx

д!3

дx

+ Л(ї\

дГ

дx

+ Ed1,

, I4) дИ=0,

дк

, I4) дИ = 0.

(2)

Умножим первое уравнение системы на

дy

а третье на —— и сложим их. В результате получим

= 0.

ду дх ду дх

Это означает, что якобиан преобразования

з=дШ=о.

д( х, у)

Следовательно, Е3 есть некоторая функция от Е1 . Аналогично из второго и третьего уравнений получаем, что Е4 есть некоторая функция от Е2.

Поэтому

ЖЕ1,..., Е4) = ЖЕ1, Е2), ЖЕ1,..., Е4) = ЖЕ1, Е2).