[ «РАЗВИТИЕ ТЕОРИИ И МЕТОДОВ ОПРЕДЕЛЕНИЯ ПРОСТРАНСТВЕННЫХ КООРДИНАТ НАЗЕМНЫХ ПУНКТОВ ПО РЕЗУЛЬТАТАМ ТРАЕКТОРНЫХ ИЗМЕРЕНИЙ СПУТНИКОВ ...»
Министерство образования и наук
и Российской Федерации
Федеральное государственное бюджетное образовательное учреждение
высшего профессионального образования
«СИБИРСКАЯ ГОСУДАРСТВЕННАЯ ГЕОДЕЗИЧЕСКАЯ АКАДЕМИЯ»
(ФГБОУ ВПО «СГГА»)
На правах рукописи
Дударев Владимир Иванович
РАЗВИТИЕ ТЕОРИИ И МЕТОДОВ ОПРЕДЕЛЕНИЯ
ПРОСТРАНСТВЕННЫХ КООРДИНАТ НАЗЕМНЫХ ПУНКТОВ ПО
РЕЗУЛЬТАТАМ ТРАЕКТОРНЫХ ИЗМЕРЕНИЙ СПУТНИКОВ
25.00.32 – “Геодезия” Диссертация на соискание учёной степени доктора технических наукНаучный консультант – доктор технических наук, профессор Карпик Александр Петрович Новосибирск –
ОГЛАВЛЕНИЕ
ВВЕДЕНИЕ.…………………………………………………………………………1 ТЕОРИЯ ОЦЕНИВАНИЯ ПАРАМЕТРОВ СОСТОЯНИЯ НЕЛИНЕЙНЫХ
ДИНАМИЧЕСКИХ СИСТЕМ….………………………………………………... 1.1 Постановка задачи оценивания состояния нелинейных динамических систем в спутниковой геодезии..……………………………………….………… 1.2 Системы координат и шкалы времени..……………………….…………….. 1.3 Преобразования координат……………………………………………………. 1.4 Общие положения и математическая постановка задачи определения геодезических параметров по результатам радиотехнических траекторных измерений спутников……….………………………………..……………………. 1.5 Выбор системы элементов орбиты спутника и фундаментальной матрицы …. 1.6 Оценка влияния изменения гравитационного поля Земли на движение космического аппарата....………………………………………………………….2 МАТЕМАТИЧЕСКИЕ МОДЕЛИ И УРАВНЕНИЯ ПОПРАВОК ДЛЯ
ОБРАБОТКИ РЕЗУЛЬТАТОВ РАДИОТЕХНИЧЕСКИХ ТРАЕКТОРНЫХ
ИЗМЕРЕНИЙ КОСМИЧЕСКИХ АППАРАТОВ………………………………... 2.1 Краткая характеристика современных методов и математических моделей радиотехнических измерений..……………………….………..………. 2.2 Влияние рефракции земной атмосферы на траекторию радиосигнала…….. 2.3 Принцип доплеровских методов определения местоположения…………… 2.4 Математические модели доплеровских траекторных измерений…………. 2.5 Математические модели радиодальномерных траекторных измерений….. 2.6 Методика учёта влияния атмосферной рефракции в радиотехнических траекторных измерениях спутников……………..…..…………….……………. 2.7 Уравнения поправок для определения геодезических параметров по результатам радиотехнических траекторных измерений спутников………….3 ВОПРОСЫ НАБЛЮДАЕМОСТИ В ЗАДАЧЕ ОПРЕДЕЛЕНИЯ
ГЕОДЕЗИЧЕСКИХ ПАРАМЕТРОВ…………………………….……….……... 3.1 Вопросы наблюдаемости в космической геодезии…………………….…… 3.2 Условия наблюдаемости параметров движения спутника при измерении его наклонной дальности и лучевой скорости………………..………..……….. 3.3 Условия наблюдаемости пространственного положения наземного пункта при измерении наклонной дальности и лучевой скорости спутника…. 3.4 Условия наблюдаемости в задаче определения параметров вращения Земли по спутниковым данным …………………………..……………..……… 3.5 Условия наблюдаемости при определении ЭВО по ГНСС-измерениям…. 3.6 Методика оценки влияния ошибок расчета матрицы коэффициентов и вектора правой части на решение СЛАУ в некоторых задачах космической геодезии…………………………………………….…….….……………………. 3.7 Точность определения высот пунктов по результатам ГНСС-измерений…4 МЕТОДЫ ПОСТРОЕНИЯ И УРАВНИВАНИЯ ГЕОДЕЗИЧЕСКИХ
СЕТЕЙ С ИСПОЛЬЗОВАНИЕМ ГНСС-ТЕХНОЛОГИЙ………..……………. 4.1 Использование ГНСС-измерений для решения геодезических задач..……. 4.2 Развитие геодезических сетей пространственными векторными построениями ………………………………………..…………………………… 4.3 Уравнивание геодезических сетей по результатам относительных ГНСС-измерений ………………………………………..……………….………. 4.4 Расчет аномалии высоты при ГНСС-определениях пространственных координат точек земной поверхности.………………………..……...…………. 4.5 Методика выявления эффекта морозного выпирания пунктов.……..……. 4.6 Некоторые приемы и особенности применения ГНСС-технологий при выполнении топографо-геодезических работ …………………..……...………. ЗАКЛЮЧЕНИЕ …………………………………………………..……..…….…. СПИСОК СОКРАЩЕНИЙ И УСЛОВНЫХ ОБОЗНАЧЕНИЙ……………….. СПИСОК ЛИТЕРАТУРЫ…………………………. ……………………..…..…. ПРИЛОЖЕНИЕ А (рекомендуемое) КЛАССИФИКАЦИЯ КООРДИНАТНЫХ СИСТЕМ…………………………….. ПРИЛОЖЕНИЕ Б (рекомендуемое) СИНГУЛЯРНОЕ РАЗЛОЖЕНИЕ ……... ПРИЛОЖЕНИЕ В (справочное) ВЫВОД АНАЛИТИЧЕСКИХ ВЫРАЖЕНИЙ ДЛЯ РАСЧЕТА ЭЛЕМЕНТОВ МАТРИЦЫ X/Y.. ПРИЛОЖЕНИЕ Г (справочное) ВЫВОД АНАЛИТИЧЕСКИХ ВЫРАЖЕНИЙ ДЛЯ РАСЧЕТА ЭЛЕМЕНТОВ МАТРИЦАНТА….. ПРИЛОЖЕНИЕ Д (рекомендуемое) ОЦЕНКА ОТНОСИТЕЛЬНОЙ ОШИБКИ МАТРИЦАНТА………..………………….....… ПРИЛОЖЕНИЕ Е (рекомендуемое) ВИД И АНАЛИЗ МОДЕЛЕЙ РАДИОТЕХНИЧЕСКИХ ТРАЕКТОРНЫХ ИЗМЕРЕНИЙ……… ПРИЛОЖЕНИЕ Ж (обязательное) ИНТЕГРАЛ КВАДРАТА СКОРОСТИ ИЗМЕНЕНИЯ ДЛИНЫ ФАЗОВОГО ПУТИ ЭМВ……… ПРИЛОЖЕНИЕ И (рекомендуемое) ЧИСЛЕННЫЕ ЗНАЧЕНИЯ ИОНОСФЕРНОЙ И ТРОПОСФЕРНОЙ ПОПРАВОК..………... ПРИЛОЖЕНИЕ К (обязательное) ВИД НЕКОТОРЫХ КОЭФФИЦИЕНТОВ В УРАВНЕНИЯХ ПОПРАВОК…...………………….…... ПРИЛОЖЕНИЕ Л (обязательное) ОШИБКИ ОПРЕДЕЛЕНИЯ ЭВО………... ПРИЛОЖЕНИЕ М (обязательное) ОШИБКИ ОПРЕДЕЛЕНИЯ ВЫСОТ НАЗЕМНЫХ ПУНКТОВ…………………………………. ПРИЛОЖЕНИЕ Н (справочное) РАЙОН САМБУРГСКОГО МЕСТОРОЖДЕНИЯ ГАЗА….……………………….…….………. ПРИЛОЖЕНИЕ П (рекомендуемое) ПРОТОКОЛ УТИЛИТЫ “Baselines" ….. ПРИЛОЖЕНИЕ Р (рекомендуемое) ПРИМЕР УРАВНИВАНИЯ ВЕКТОРНОГО ХОДА…………………...………………… ПРИЛОЖЕНИЕ С (рекомендуемое) ПРИМЕР УРАВНИВАНИЯ ВЕКТОРНОЙ ЗАСЕЧКИ………………….…..…………... ПРИЛОЖЕНИЕ Т (рекомендуемое) ПРИМЕР ФОРМИРОВАНИЯ СИСТЕМЫ ЛИНЕЙНЫХ УРАВНЕНИЙ ПОПРАВОК….………ВВЕДЕНИЕ
Актуальность темы исследования. Растущий интерес человечества к вопросам детального изучения физической поверхности Земли, геодинамических процессов, совершенствование аэрокосмических методов и средств морской, воздушной и наземной навигации, а также необходимость решения целого ряда научных и практических задач потребовали создания новых технологий выполнения геодезических работ и методов оперативной математической обработки результатов измерений, полученных с помощью навигационных спутников и их систем. Высокая точность решения задач космической геодезии возможна только при развитии теории и методов математической обработки больших массивов результатов траекторных измерений спутников для определения необходимых орбитальных и геодезических параметров. В данном исследовании особое внимание уделяется: параметрам движения спутника, параметрам вращения Земли (координатам истинных полюсов Земли и поправке к неравномерной шкале Всемирного времени UT1), элементам взаимного ориентирования (ЭВО) геодезических систем координат, пространственным координатам наземных пунктов (НП), под которыми понимаются точки различного назначения, расположенные как на земной и водной поверхности, так и в воздушном пространстве.Современная спутниковая радиотехническая аппаратура позволяет измерять длину траектории распространения электромагнитного колебания от космического аппарата (КА) до НП с точностью до 0,01длины волны. Это дает возможность гарантированно определять орбитальные и геодезические параметры по результатам радиотехнических траекторных измерений спутников с высокой точностью при различной геометрии относительного расположения КА и НП.
Чтобы решать с такой высокой точностью поставленную задачу, требуются адекватные этим измерениям математические модели и соответствующие им уравнения поправок. Для этого требуется развивать теорию и методы, связанные с: учетом влияния на результаты измерений атмосферы Земли, особенностей работы приемной и передающей радиотехнической спутниковой аппаратуры и других физических эффектов; разработкой оптимальных математических алгоритмов формирования матриц коэффициентов и вектора правой части уравнений поправок; оценкой влияния изменения гравитационного поля Земли на движение КА.
Существующие математические методы уравнивания плановых и высотных геодезических сетей, созданных посредством выполненных относительных ГНССизмерений, используют не все потенциальные возможности геодезических навигационных спутниковых систем типа ГЛОНАСС, НАВСТАР и др. Поэтому развитие теории и методов определения координат одного или нескольких НП при использовании различных пространственных геодезических построений на основе относительных ГНСС-измерений имеет важное научное значение.
Задача преобразования высот точек из геодезической системы в нормальную и обратно имеет большое научное и практическое значение. Существующие математические методы решения этой задачи достаточно точны, но отличаются своей громоздкостью и большим объемом подготовительных и вычислительных работ: требуют построения так называемого “локального геоида” или интерполирования аномалии высоты между исходными НП. Кроме того, эти методы, как правило, ориентированы на значительные по площади территории.
При выполнении геодезических работ с использованием ГНСС-технологий на небольших территориях требуется оперативное решение такой задачи. Поэтому разработка теории оптимального преобразования высот точек из геодезической системы в нормальную и обратно является актуальной задачей.
Степень разработанности темы. За предыдущие годы накоплен достаточно большой опыт и объем знаний в области определения орбитальных и геодезических параметров по результатам траекторных измерений спутников.
Весомый вклад в развитие этого направления внесен научными сотрудниками Центрального научно-исследовательского института геодезии, аэрофотосъемки и картографии (ЦНИИГАиК); Научно-исследовательского института Военнотопографической службы (НИИ ВТС); Московского государственного университета геодезии и картографии (МИИГАиК); Открытого акционерного общества «Информационные спутниковые системы» им. академика М. Ф. Решетнева»;
Сибирской государственной геодезической академии (НИИГАиК), Института астрономии РАН (Астрономический Совет АН СССР); Государственного астрономического института им. П. К. Штенберга и целого ряда других организаций. Значительный вклад в решение научных и практических задач космической геодезии привнесли советские, российские и зарубежные ученые, в том числе: М. Бурша, Г. Вейс, Г. Н. Дубошин, Б. Ф. Жданюк, А. А. Изотов, М. С. Урмаев, П. Е. Эльясберг, У. Каула и др.
Для обработки результатов наблюдений спутников отечественными и зарубежными исследователями разработан ряд математических моделей радиотехнических траекторных измерений. При анализе этих моделей было выявлено их несоответствие возросшему за последние годы уровню точности измерений, выполняемых современной спутниковой аппаратурой. По этой причине возникла потребность в создании более совершенных математических моделей, учитывающих влияние на результаты измерений как внешних, так и внутренних факторов, обусловленных атмосферой Земли, спецификой работы приемной и передающей спутниковой радиотехнической аппаратуры, особенностями движения КА по орбите и др. Разработка более точных моделей траекторных измерений обусловила и разработку соответствующих им уравнений поправок, которые должны обеспечивать надежную оценку определяемых параметров.
Теоретическим разработкам в области математической постановки решения задачи оценивания посвящены работы В. Н. Брандина, Б. Ф. Жданюка, Г. Н. Разоренова, И. В. Онькова, П. Е. Эльясберга и ряда других авторов.
Повышение точности и увеличение видов траекторных измерений спутников, рост количества определяемых параметров требует постоянного совершенствования математических алгоритмов обработки результатов измерений, что вызывает необходимость в дальнейшем развитии и совершенствовании математической постановки решения задачи оценивания.
В вопросах развития и уравнивания геодезических сетей с использованием ГНСС-технологий накоплен значительный опыт. Решению таких задач посвящены Ю. И. Маркузе, Е. Г. Бойко, А. В. Антипова, Г. А. Шанурова и др. Особого внимания заслуживает работа Ю. И. Маркузе и W. M. Welsh, в которой представлена система линейных уравнений поправок для уравнивания геодезических сетей по результатам относительных ГНСС-измерений. Для обработки ГНСС-измерений используются программные комплексы, например:
“GPSurvey”, “Trimble Geomatics Office”, “Ski” и др. Однако, в руководствах по их применению не приводится описание математических алгоритмов обработки этих измерений, что существенно усложняет осмысление получаемых с их помощью результатов. Кроме этого, подобные зарубежные программные продукты не позволяют активно влиять на ход вычислительного процесса. Универсальный и оптимальный математический алгоритм решения рассматриваемой задачи, повидимому, создать невозможно, так как он будет слабо учитывать особенности отдельно поставленной задачи. Поэтому существует потребность в разработке математических алгоритмов, позволяющих эффективно решать конкретные задачи с учетом их специфики. К ним относятся задачи по определению координат как одного или двух, так и нескольких наземных пунктов в системе координат исходных геодезических пунктов, элементов взаимного ориентирования, преобразованию высот точек из геодезической системы в нормальную и обратно, по созданию геодезических сетей – сплошных, протяженных, высотных и других.
Цели и задачи исследования. Целью исследования является развитие теории и методов гарантированного определения по результатам современных беззапросных радиотехнических траекторных измерений спутников пространственных координат наземных пунктов и орбитальных параметров с точностью, соответствующей точности современных спутниковых радиотехнических систем (РТС).
Для достижения поставленной цели необходимо было решить следующие основные задачи:
– систематизировать и развить теоретическое обоснование регулярной постановки задачи оценивания геодезических параметров по результатам беззапросных радиотехнических траекторных измерений КА;
– разработать теорию формирования нелинейных математических моделей беззапросных радиотехнических траекторных измерений спутников, отвечающих современному уровню точности;
– разработать уравнения поправок для обработки результатов беззапросных радиотехнических траекторных измерений спутников;
– развить теорию и методику учета влияния атмосферной рефракции на результаты радиотехнических траекторных измерений спутников;
– разработать аналитическую теорию оценки влияния изменения гравитационного поля Земли (ГПЗ) на движение спутника по орбите;
– развить теоретические положения и выполнить экспериментальные исследования для выявления такой геометрии взаимного расположения КА и НП, при которой по выполненным радиотехническим траекторным измерениям определяются не все оцениваемые параметры;
– развить теорию уравнивания и методы построения различных геодезических сетей по результатам относительных ГНСС-измерений;
– разработать оптимальный математический алгоритм преобразования высот точек из геодезической системы в нормальную и обратно.
Объект исследования – динамическая система, состоящая из нескольких космических аппаратов и наземных пунктов, между которыми выполняются беззапросные радиотехнические траекторные измерения.
Предмет исследования – теория и методы определения пространственных координат наземных пунктов и параметров движения спутников по результатам траекторным измерений.
Научная новизна:
– разработанные нелинейные математические модели и уравнения поправок для обработки результатов беззапросных радиотехнических траекторных измерений спутников позволяют получать оценки заданного набора определяемых орбитальных и геодезических параметров при различной геометрии взаимного расположения КА и НП с точностью, соответствующей точности выполняемых измерений;
– предложенные методы построения и уравнивания различных геодезических сетей по результатам относительных ГНСС-измерений уменьшают время и производственные затраты на выполнение топографо-геодезических работ;
– разработанный математический алгоритм преобразования высот точек из геодезической системы в нормальную и обратно повышает эффективность выполнения топографо-геодезических работ.
Теоретическая и практическая значимость работы. Разработанные в диссертации теоретические и методические положения позволяют оптимизировать решение измерительной задачи, связанной с определением орбитальных параметров и пространственных координат геодезических пунктов.
Полученные аналитические выражения доведены до практического использования и предназначены для решения фундаментальных и прикладных задач геодезии, космической геодезии, геодинамики и точной навигации.
Результаты исследований предлагается использовать для координатновременного обеспечения целого ряда задач народно-хозяйственного и оборонного назначения. К ним относятся: изучение топографии земной поверхности и геодинамических процессов, поиск и освоение полезных ископаемых, геодезическое обеспечение строительства различных инженерных сооружений, определение параметров движения КА. Результаты исследований позволяют повысить эффективность проведения производственных работ и существенно снизить затраты на их выполнение.
Методология и методы исследования основаны на системном анализе, теории и математических методах изучения поведения нелинейных динамических систем, теоретических исследованиях в области космической геодезии и баллистики, небесной механики и астродинамики, астрометрии, теории фигуры и физики атмосферы Земли. Решение поставленных задач базировалось на использовании математического моделирования, анализа и синтеза, сравнения, обобщений и оценок с привлечением методов интегрального и дифференциального исчисления, численных методов линейной алгебры, теории математической обработки и интерпретации результатов геодезических измерений.
Положения, выносимые на защиту:
– результаты теоретических разработок в области формирования нелинейных и линеаризованных математических моделей радиотехнических траекторных измерений спутников, отвечающих уровню точности современных РТС, повышают степень состоятельности оценок определяемых параметров;
– теория расчета различия между длиной траектории распространения радиосигнала в тропосфере Земли и прямой наклонной дальностью до спутника повышает степень адекватности математических моделей радиотехнических траекторных измерений спутников;
– аналитическая теория учета влияния изменения гравитационного поля Земли на движение спутника по орбите обеспечивает существенное преимущество перед численным методом по эффективности, быстродействию машинных алгоритмов и надежности получаемых результатов;
– теоретические и экспериментальные исследования условий наблюдаемости орбитальных и геодезических параметров по результатам радиотехнических траекторных измерений спутника повышают степень регулярной постановки измерительной задачи и способствуют созданию оптимальных алгоритмов ее решения;
– теоретические разработки и методика оптимального преобразования высот точек из геодезической системы в нормальную и обратно, а также развитая теория уравнивания и методы развития различных геодезических сетей по результатам относительных ГНСС-измерений, существенно повышают эффективность использования ГНСС-измерений и понижают производственные затраты при проведении геодезических работ.
Соответствие диссертации паспорту научной специальности. Тематика диссертации соответствует следующим пунктам областей исследований паспорта научной специальности 25.00.32 «Геодезия»: 2 – «Создание геодезической координатно-временной основы различного назначения с использованием геодезических, астрономических, гравиметрических и других (космических, наземных и подземных) методов измерений; оценка их степени устойчивости и характера изменений, вопросы их проектирования и оптимизации. Геодезические системы координат»; 4 – «Разработка новых принципов, методов, технических средств и технологий геодезических измерений для определения геометрических и физических параметров Земли, ее поверхности, объектов, явлений и процессов на ней, в том числе для производства наземных топографических съемок»;
11 – «Теория и практика математической обработки результатов геодезических измерений и информационное обеспечение геодезических работ. Автоматизированные технологии создания цифровых трехмерных моделей технологических объектов, процессов и явлений по геодезическим данным».
Степень достоверности и апробация результатов. Основные положения и результаты диссертационной работы докладывались, обсуждались и получили одобрение на: Всесоюзном научно-организационном совещании ИТА АН СССР (г. Ленинград, 1990 г.); IX Съезде ВАГО (г. Новосибирск, 1990 г.); 42 и 43 научнотехнических конференциях преподавателей НИИГАиК (г. Новосибирск, 1993, 1994 г.г.); 44 научно-технической конференции с международным участием (г. Новосибирск, 1994 г.); 45 – 48 научно-технических конференциях преподавателей СГГА (г. Новосибирск, 1995 – 1998 гг.); Третьем Сибирском конгрессе по прикладной и индустриальной математике (ИНПРИМ-98) (г. Новосибирск, 1998 г.); Международной научно-технической конференции, посвященной 220-летию МИИГАиК (г. Москва, 1999 г.); 50 научно-технической конференции преподавателей СГГА (г. Новосибирск, 2000 г.); Четвертом Сибирском конгрессе по прикладной и индустриальной математике (ИНПРИМ-2000) (г. Новосибирск, 2000 г.); 51 научно-технической конференции преподавателей СГГА, посвященной памяти академика В. В. Бузука (г. Новосибирск, 2001 г.); 53 Международной научно-технической конференции, посвященной 70-летию СГГА (г. Новосибирск, 2003 г.); 54 научно-технической конференции, посвященной 225-летию геодезического образования в России (г. Новосибирск, 2004 г.); Международной (г. Москва, 2004 г.); 2-й Региональной научно-практической конференции (г. Иркутск, 2006 г.); Международных научных конгрессах “ГЕО-СИБИРЬ” (г. Новосибирск, 2006 – 2008 гг., 2010 – 2012 гг.).
Работа выполнялась в рамках: госбюджетной научно-исследовательской работы (НИР) “Анализ и разработка алгоритмов оценивания геодезических и геодинамических параметров по результатам траекторных наблюдений ИСЗ” (1992 г. – 1996 г.), фундаментальной НИР № ГР 0186.0036274 “Разработка методов изучения динамики гравитационного поля и фигуры Земли планетарного, регионального и локального характера” (1995 г.) по комплексной научнотехнической программе “Человек и окружающая среда”, фундаментальной НИР № ГР 0196.0005075 Грант Госкомвуза РФ “Исследование динамической составляющей гравитационного потенциала Земли и его характеристик, построение их математических моделей” (1996 г.), хоздоговорной НИР № 4.30.392 “Обследование приаэродромной территории, определение координат и высот препятствий” (2001 г.), госбюджетной НИР № 6.30.194 “Развитие локальных и региональных геодезических сетей с использованием результатов спутниковых траекторных измерений” (2002 г. – 2006 г.), фундаментальной НИР № ГР 012008.03158 “Разработка теории и методов определения геодезических и геодинамических параметров по результатам спутниковых радиотехнических траекторных измерений” (2008 г. – 2010 г.) по аналитической ведомственной целевой программе “Развитие вузовского потенциала высшей школы России на 2009 – 2012 гг.”, фундаментальной НИР № ГР 01201151296 “Исследование и разработка теории и методов оценки состояния нелинейных динамических систем по результатам траекторных измерений космических аппаратов” (2011 г.) по аналитической ведомственной целевой программе “Развитие вузовского потенциала высшей школы России на 2009 – 2012 гг.”, фундаментальной НИР № ГР 01201252771 государственного задания Министерства образования и науки РФ “Теория и методы обработки радиотехнических траекторных измерений космических аппаратов” (2012 г.).
Разработки автора нашли применение при выполнении производственных работ, связанных с развитием геодезических сетей и проведением топографических съемок участков земной поверхности с использованием ГНССтехнологий, в ряде регионов России.
Результаты выполненных исследований внедрены в экспедицию № Ноябрьского филиала ФГУП “ЗапсибАГП”, в ФГУП Кузбасское топографомаркшейдерское предприятие “КУЗБАССМАРКШЕЙДЕРИЯ”, в учебный процесс СГГА.
Публикации по теме диссертации. Основные результаты диссертационных исследований представлены в 36 публикациях, из которых 12 – в журналах, входящих в перечень российских рецензируемых научных журналов и изданий.
Структура диссертации. Диссертация состоит из введения, четырех разделов, заключения, списка литературы из 195 наименований. Общий объем диссертации составляет 314 страниц машинописного текста, содержит 22 рисунка, 5 таблиц, 16 приложений.
1 ТЕОРИЯ ОЦЕНИВАНИЯ ПАРАМЕТРОВ СОСТОЯНИЯ НЕЛИНЕЙНЫХ
ДИНАМИЧЕСКИХ СИСТЕМ
1.1 Постановка задачи оценивания состояния нелинейных динамических систем в спутниковой геодезии Задачи спутниковой геодезии относятся к классу измерительных задач, в которых по результатам измерений требуется определить состояние конкретной динамической системы на заданном отрезке времени T. Такие задачи иначе называются обратными задачами или задачами оценивания [125, 141, 161].Наличие больших массивов измерительной информации, а также желание необходимость создания оптимальных алгоритмов (в смысле их универсальности, регулярности и корректности, а также оперативности решения задачи) программного обеспечения систем автоматизированной математической обработки и интерпретации результатов измерений.
В спутниковой геодезии измерительная информация получается в ходе выполнения траекторных наблюдений КА с помощью специальной аппаратуры в течение одного или нескольких сеансов. Их продолжительность и промежутки времени между ними, а также распределение измерений внутри сеансов зависит от условий проведения измерений. В ряде случаев их повторить не удается. Поэтому планирование проведения наблюдений и постановка решения задачи оценивания требуют к себе особого внимания и подхода.
математической модели динамической системы (ДС), уравнений измерений, условий опыта, критерия качества и метода решения нелинейной системы уравнений [14, 88, 162]. Рассмотрим общие вопросы постановки задачи оценивания состояния нелинейной ДС, состоящей из нескольких КА и НП, по результатам беззапросных радиотехнических траекторных измерений спутников [38, 85, 133].
Под динамической системой понимается множество элементов, изменяющих свое поведение во времени и объединенных между собой в единое целое совокупностью некоторых соответствий и отношений [14, 15, 88]. Состояние реальной системы полностью описывается бесконечным числом параметров. В условиях обратной задачи имеется конечное число измерений, искаженных ошибками.
Поэтому в большинстве случаев фактическое состояние этой системы определить не удается. В связи с этим она заменяется некоторой математической моделью, зависящей от конечного числа параметров. Их совокупность в текущий момент времени tT образует r-мерный вектор-столбец параметров состояния q = q(t):
динамической системы при решении научных и практических задач. Модели должны быть наиболее полными, чтобы однозначно и достаточно точно описывать поведение реальной системы, относительно простыми и удобными для вычислений. Задача оценивания подразумевает определение некоторой совокупности параметров q ДС, характеризующих ее состояние в некоторый заданный момент времени tT. В спутниковой геодезии под ним понимается так называемый начальный момент времени t0T. Часто t0 находится примерно посередине временного отрезка Т.
Начальное состояние ДС в момент t0 определяется вектором q0 = q(t0) начальных условий (НУ). Для прогнозируемых моделей текущее состояние системы устанавливается либо соотношением либо после интегрирования обыкновенных дифференциальных уравнений В общем случае вектор-функции F и f нелинейные. Зависимости (1.2) и (1.3) выражают основные закономерности процесса движения (динамики поведения системы) и называются моделями движения [14, 88]. Они подразделяются на геометрическое описание движения, в динамических моделях для описания процесса движения учитываются сведения о массе и воздействующих силах [160].
Непосредственно измеряемую величину Z i' (разность фаз, частота и т.д.) принято называть измеряемым параметром. Так как в один момент времени t могут выполняться несколько видов радиотехнических траекторных измерений КА (интегральные и дифференциальные доплеровские, радиодальномерные и другие наблюдения), то индекс i будет принимать значения i = 1, 2,..., d (d – число одновременно измеряемых параметров). Поэтому эти параметры в момент t образует d-мерный вектор-столбец Z' = Z' (t) измеряемых параметров:
В результатах измерений могут присутствовать систематические ошибки.
Будем рассматривать такие систематические ошибки, значения которых меняются не только между сеансами измерений, но и в пределах одного отдельно взятого сеанса. Такого рода ошибки называются сингулярными [13, 15]. Их появление обусловлено воздействием на результаты измерений ряда медленно меняющихся со временем факторов. К ним относятся аберрационные и релятивистские эффекты, рефракция атмосферы Земли, а также мешающие параметры. Эти параметры практически постоянны в пределах одного сеанса измерений и их появление обусловлено нестабильностью работы генераторов передающей и приемной радиотехнической аппаратуры, внутренними временными задержками передатчика и приемника, нестабильностью состояния атмосферы Земли в течение одного сеанса радиотехнических измерений. Этот вопрос будет подробно рассмотрен в разделе 2. Для компонент вектора мешающих параметров Г = Г(t) [162] отсутствуют математические модели. В связи с вышесказанным под измеряемыми псевдопараметрами Z будем понимать измеряемые параметры Z', беззапросных радиотехнических траекторных измерений под измеряемыми псевдопараметрами будем понимать псевдодальность, разность псевдодальностей и псевдоскорость.
Функциональная зависимость определяемая в общем случае нелинейной измеряемой d-мерной вектор-функцией на один текущий момент времени t, представляет собой матричную запись уравнений измерений и является математической моделью измерений по их видам (импульсные и фазовые радиодальномерные, дифференциальные и интегральные доплеровские и т. д.) [13 – 15, 88] В ней параметр является k-мерным вектор-столбцом мешающих параметров для одного отдельно взятого сеанса измерений. Так как сеансы достаточно короткие, то будем считать вектор Г постоянным в течение одного сеанса, но изменяющимся случайным образом от сеанса к сеансу.
Результаты измерений Z i (i = 1, 2,..., d; d – число одновременно выполненных видов измерений) на момент времени t, формируют d-мерный вектор-столбец и содержат случайные ошибки измерений i (i = 1, 2,..., d), образующие случайный d-мерный вектор-столбец = (t) Вектор измерений Z функционально связан с d-мерным вектором-столбцом истинных значений измеряемых параметров Z = [Z1 Z 2...Z d ]T и случайным вектором ошибок измерений. Эта зависимость описывается нелинейной d-мерной вектор-функцией в которой параметров состояния q = q(t ), определяющего поведение реальной системы в текущий момент времени t. Способ комбинации ошибок измерений с измеряемыми параметрами и статистические свойства этих ошибок в математической статистике называются условиями опыта. Будем рассматривать способ комбинации ошибок измерений с вектором Z, при котором (Z, 0, t ) = Z.
Такому условию удовлетворяет аддитивный способ комбинации ошибок [14] Совместно со случайными ошибками измерений на точность решения задачи оценивания оказывают влияние и ошибки модели измерений (1.5). Их появление обусловлено ограниченным числом принятых параметров состояния ДС. Эти ошибки в момент времени t образуют случайный d-мерный вектор-столбец =(t):
С учетом этих ошибок вектор-функция в выражении (1.11) примет вид:
Тогда вектор-столбец Z результатов измерений в равенстве (1.12) можно определить как То есть результаты измерений могут быть представлены соответствующими математическими моделями измерений, искаженных случайными ошибками измерений и моделей измерений.
времени [t1,t2,…, th]T (h – число измерений в сеансе), совокупные моменты векторы-столбцы Z первичных измеряемых псевдопараметров и Z результатов измерений размерности dxh определим следующими зависимостями:
В них Здесь вектор-столбец Z образует выборку результатов измерений одного сеанса.
Наиболее полно векторы и представляются многомерными законами распределения, которые практически невозможно записать. Поэтому на практике обычно задаются первыми двумя моментами многомерного распределения ошибок: математическими ожиданиями M( ), M() и ковариационными матрицами D( ), D(). Для нормального закона распределения ошибок функция плотности распределения вероятностей полностью определяется этими характеристиками.
Как уже отмечалось выше, помимо случайных ошибок и на точность решения измерительной задачи оказывают влияние и сингулярные ошибки. Они функционально связаны с вектором мешающих параметров Г, который обусловлен специфическими условиями выполнения радиотехнических траекторных измерений КА.
Разложив правую часть выражения (1.5) в ряд Тейлора, ограничившись линейными членами, запишем Здесь частная производная есть прямоугольная dxk матрица частных производных от измеряемой векторфункции вида (1.6) по компонентам вектора мешающих параметров, вычисленная при значении Г=0. Выражение (1.19) представляет собой линеаризованные в окрестности векторов Г=0 и q уравнения измерений (1.5). В них сингулярные ошибки =(t) по видам измерений, образующие d-мерный вектор-столбец рассчитываются по формуле На основании зависимостей (1.19) и (1.22) уравнения измерений (1.5) с аддитивным способом комбинации вектора-столбца Z = (q, =0, t ) первичных измеряемых параметров и сингулярных ошибок на момент времени t могут быть записаны как Выражение (1.23) показывает связь вектора Z измеряемых псевдопараметров и вектора Z' измеряемых параметров.
Матрица П в равенстве (1.20) зависит от времени. Следовательно, сингулярные ошибки также зависят от времени и изменяются не только между сеансами измерений, но и в пределах одного сеанса. Влияние этих ошибок на результаты измерений выборки, сформированной из одного сеанса, представим выражением в котором В равенстве (1.24) сингулярные ошибки образуют совокупный вектор-столбец размерности d·h вида Вектор-столбец для d-видов радиотехнических измерений в одном сеансе можно определить по формуле в которой матрица есть прямоугольная матрица размерности (d·h) x k вида В спутниковой геодезии под вектором Z обычно понимают так называемые радионавигационного сигнала, разность фаз принятого и опорного колебаний, время распространения радиоимпульса от КА до НП и т.д. В большинстве случаев математическую обработку результатов измерений удобнее выполнять с геометрический и физический смысл – скорость относительного перемещения КА и НП, наклонная дальность от НП до КА и т.д. Под вторичными измеряемыми псевдопараметрами для радиотехнических траекторных измерений КА псевдоскорость относительного перемещения КА и НП, псевдодальность между НП и КА и т.д. Иначе говоря, за вторичные измеряемые псевдопараметры принимаются вторичные измеряемые параметры, искаженные сингулярными ошибками [74]. Вторичные измеряемые псевдопараметры в момент времени t также образуют d-мерный вектор-столбец вида Между первичными Z и вторичными измеряемыми псевдопараметрами существует функциональная зависимость в которой W в общем случае является нелинейной первичной измеряемой d-мерной вектор-функцией, вид которой определяется типом измерений. В свою очередь, на момент времени t вектор определяется нелинейной вторичной измеряемой d-мерной вектор-функцией от текущих значений вектора q параметров состояния и вектора Г мешающих параметров:
Разложив вектор-функцию из равенства (1.30) в ряд Тейлора в окрестности векторов Г=0 и q, ограничившись при этом линейными членами разложения, запишем аналитическую зависимость между вторичными измеряемыми псевдопараметрами и их сингулярными ошибками Здесь частная производная есть прямоугольная dxk матрица частных производных от вектор-функции по компонентам вектора мешающих параметров, вычисленная при значении Г=0. В этом случае систему уравнений (1.29) можно записать в виде, в котором вторичные измеряемые параметры '=(q,Г=0,t) и их сингулярные ошибки имеют аддитивный способ комбинации:
В этой системе уравнений ( ) есть d-мерная матрица вторичных измеряемых псевдопараметров. Вектор-столбец сингулярных ошибок вторичных измеряемых параметров рассчитывается по формуле Матрица Е зависит от времени. Следовательно, сингулярные ошибки также являются функциями времени и изменяются в пределах одного сеанса измерений.
Зависимость между результатами измерений, выполненных в течение одного сеанса, и их сингулярными ошибками имеет вид в которой определяется из произведения прямоугольной (d·h)xk матрицы E и k-мерного вектора-столбца мешающих параметров Г:
Для моделей измерений, представленных зависимостями (1.23) и (1.33), матрицы частных производных П и E могут быть найдены как При обработке измерений, полученных в течение нескольких сеансов, формируются уравнения измерений вида (1.24) и (1.36) по всем сеансам.
Целью измерительной задачи является получение несмещенных оценок определяемых параметров. Для выполнения условия несмещенности требуется, чтобы математические ожидания ошибок траекторных измерений, моделей измерений и сингулярных ошибок (по всем сеансам) равнялись нулю, то есть M( ) = 0, M() = 0, M() = 0 (либо M() = 0 ). Так как то для выполнения требования M() = 0 (или M() = 0 ) должно выполняться равенство В общем случае математическое ожидание вектора мешающих параметров Г не равняется нулю, поэтому его компоненты следует включать в состав вектора оцениваемых параметров [13, 15]. Это позволяет повысить точность найденных значений определяемых параметров не менее, чем на 10 % [119]. При этом не требуется каких-либо априорных сведений о вероятностных характеристиках вектора Г [162].
Выбор критерия качества на этапе постановки решения задачи оценивания определяет оптимальность получаемых оценок q' вектора параметров состояния нелинейной динамической системы. Это может быть метод наименьших квадратов, метод наименьших модулей и т. д. Метод статистической обработки результатов измерений обычно получает свое название по наименованию используемого критерия качества [13 – 15, 88, 162].
Математическая постановка решения измерительной задачи должна обеспечивать единственное решение с требуемыми предельными свойствами по объему выборки измерений. Для этого необходимо, чтобы ее постановка была регулярной.
Математическая постановка задачи будет регулярной, если выполняются условия адекватности, наблюдаемости и состоятельности оценок [13 – 15]. Под адекватностью понимается близость математической модели к поведению реальной ДС, под наблюдаемостью – наличие взаимно однозначного соответствия между параметрами состояния системы и ее измеряемым выходом. Состоятельность оценок подразумевает сходимость по вероятности P получаемых в процессе решения измерительной задачи оценок q'0 параметров состояния ДС к их истинным значениям q0 при сколь угодно малой > 0 и числе измерений n, стремящемся к бесконечности:
Нерегулярная постановка задачи оценивания приводит либо к множественности, либо к ошибочности ее решения.
Математическая постановка решения измерительной задачи должна быть также корректной и оптимальной. Постановка решения измерительной задачи корректная, если в процессе ее выполнения обеспечивается устойчивость решения относительно исходных данных [15, 141]. Постановка измерительной задачи будет оптимальной, если среди множества других возможных постановок ей соответствует наиболее простой и эффективный алгоритм решения. Оптимальность подразумевает также соблюдение условий регулярности и корректности.
С позиций системного подхода измерительная задача может быть представлена некоторой сложной решающей системой, у которой входом является выборка результатов измерений Z, а выходом – оценки компонент вектора q параметров состояния на некоторый начальный момент времени t0T [15]. Такая система имеет свои специфические свойства. К основным из них относятся: наблюдаемость, декомпозируемость, устойчивость, адекватность, состоятельность.
Системный подход к задаче оценивания заключается в выявлении ее характерных свойств. Затем система делится на части (подсистемы) таким образом, чтобы эти свойства сохранились для каждой ее части. После этого самостоятельному исследованию подвергается каждая подсистема в отдельности со своим входом и выходом. Информация, полученная после анализа каждой подсистемы, обобщается с целью накопления знаний о системе в целом. Кратко опишем основные свойства решающей системы, принимая во внимание тот факт, что свойства наблюдаемости, адекватности и состоятельности определяют понятие регулярности измерительной задачи, о чем было сказано выше.
Способность системы к делению на независимые части (подсистемы) называется декомпозируемостью. Вопросы декомпозиции достаточно подробно освещены в специальной литературе, в частности [14, 15].
Под устойчивостью понимается такое свойство системы, при котором малые изменения начальных условий приводят к малым изменениям измеряемого выхода [141, 142]. Следует иметь в виду, что свойство устойчивости соответствует понятию корректности постановки решения измерительной задачи.
Особо заметим, что все приведенные в данном подразделе положения, понятия и определения будут использованы при изложении последующего материала, посвященного обработке спутниковых радиотехнических траекторных измерений для определения некоторых геодезических и геодинамических параметров.
1.2 Системы координат и шкалы времени При решении геодезических задач методами космической геодезии применяются разнообразные системы координат [1, 8, 75, 89, 135], в которых математический алгоритм решения поставленной задачи оценивания будет наиболее простым. Все системы координат можно различить между собой по виду, расположению начала отсчета, выбору основной координатной плоскости, ориентировке оси абсцисс, типу.
Согласно приведенной в приложении А классификации любая выбранная система координат для полного ее определения в своем названии будет содержать все классификационные признаки [75]. Например, система координат Красовского 1942 года будет называться так: средняя земная эллипсоидальная экваториальная квазигеоцентрическая система координат. Ниже рассмотрим только те системы координат и их преобразования, которые будут применяться в настоящей работе.
Пространственное положение НП удобнее всего определять в системах координат, жестко связанных с вращающейся Землей. Это так называемые земные системы координат. В ряде случаев в них бывает также удобно определять пространственные положения естественных и искусственных небесных тел. У этих систем ось аппликат совмещается с осью вращения Земли и направляется в сторону ее северного полюса. В качестве основной координатной плоскости принимается плоскость земного экватора, а за начальный меридиан – меридиан Гринвича. Наиболее часто используются прямоугольные экваториальные геоцентрические координаты, связанные с ОЗЭ и РЭ – это общеземная (гринвичская) и референцная (геодезическая) системы.
Начало средней земной (общеземной или гринвичской) системы координат расположено в центре масс Земли, ось аппликат ZG направлена к (ОXYZ)G среднему полюсу P Земли 1900 – 1905 гг. Ось абсцисс XG направлена в точку пересечения среднего гринвичского меридиана со средним земным экватором 1900 – 1905 гг. [99, 136, 144]. За средний меридиан Гринвича принимается плоскость, содержащая вектор силы тяжести в обсерватории Гринвича и параллельная средней оси вращения Земли. Плоскость среднего меридиана Гринвича не содержит среднюю ось вращения Земли.
Начало референцной системы координат (ОXYZ)Г совмещается с центром выбранного РЭ. Ось аппликат ZГ совпадает с осью вращения этого эллипсоида. За ось абсцисс XГ принимается линия пересечения плоскости экватора РЭ и плоскости начального гринвичского меридиана, заданного геодезическими датами в исходном пункте. Ось ординат YГ дополняет систему до правой [89].
прямоугольные горизонтные системы координат. Одна из них – это средняя прямоугольная геоцентрическая горизонтная система координат (ОXYZ)H. Ее начало расположено в центре масс Земли. Ось ординат YH совпадает с линией пересечения плоскости местного меридиана НП и плоскости, проходящей через центр масс Земли параллельно плоскости местного горизонта. Ось абсцисс XH лежит в плоскости среднего экватора и дополняет систему до правой.
В описанных выше отсчетных системах положение точки остается неизменным относительно поверхности Земли и не изменяется с течением времени (если не принимать во внимание геодинамические явления).
Истинная ось вращения в теле Земли не сохраняет своего пространствен-ного положения, а испытывает периодические колебания. Следовательно, оба полюса Земли, являющиеся точками пересечения оси вращения с физической поверхностью Земли, также меняют свое положение. Отсюда следует, что земные координаты являются функцией времени. Явление перемещения мгновенных полюсов Земли относительно ее физической поверхности называется движением полюсов. Амплитуда этих перемещений доходит до 13 метров [136].
Международная служба вращения Земли (МСВЗ, IERS) отслеживает эти движения и публикует координаты xp, yp истинного полюса относительно его среднего положения P, определенного за период 1900 – 1905 гг. В связи с этим различают истинную и среднюю земные системы координат.
Начало истинной земной системы координат (ОXYZ)'G расположено в центре масс Земли. Ось аппликат Z'G направлена к истинному полюсу P' Земли в некоторый текущий момент времени tT. Ось абсцисс X'G направлена в точку пересечения истинного меридиана Гринвича с истинным экватором на эпоху t. За истинный (мгновенный) меридиан Гринвича принимается плоскость, содержащая истинную (мгновенную) ось вращения Земли и линию пересечения среднего экватора со средним меридианом Гринвича.
Обсерваториями России (Государственной системы определения параметров вращения Земли – ГС ПВЗ), Украины, Узбекистана, Польши, Болгарии и Чехии выполняются регулярные астрооптические и спутниковые наблюдения с целью определения ПВЗ: координат мгновенного полюса Земли xp, yp и Всемирного времени UT1. Окончательные данные о Всемирном времени и координатах полюса, вычисленные в Главном метрологическом центре Государственной службы времени и частоты (ГМЦ ГСВЧ), периодически публикуются в бюллетене “Всемирное время и координаты полюса” серии “Е”. Там же приводятся сведения об отличии этих данных от результатов, полученных МСВЗ. В бюллетене "Всемирное время и координаты полюса" серии "А" публикуются предварительные (срочные) значения ПВЗ и их прогноз на семь недель, а также окончательные значения ПВЗ на прошедшие даты.
Координаты xp и yp истинного полюса P' Земли определяются в касательной плоскости (координатной плоскости), проведенной к среднему полюсу P Земли.
Начало координат совмещено со средним полюсом P. Ось абсцисс xp совпадает с линией пересечения этой координатной плоскости с плоскостью среднего меридиана Гринвича. Ось ординат yp дополняет систему до правой (направлена к западу относительно меридиана Гринвича).
Для описания пространственного движения КА и НП относительно центра масс Земли наиболее приемлемой является инерциальная система координат, не вращающаяся вместе с Землей. В инерциальной системе координат начало помещается в некоторой точке пространства, либо перемещается с постоянной скоростью, направление осей в пространстве при этом сохраняется неизменным. В качестве таковой примем прямоугольную геоцентрическую экваториальную звездную систему координат (oxyz)', начало которой совмещено с центром масс Земли [136, 140]. Ось аппликат z' совпадает с ее истинной осью вращения. Ось абсцисс x' направлена в истинную точку весеннего равноденствия, которая лежит в плоскости истинного экватора в эпоху t0. Обычно за эту эпоху принимается момент задания начальных условий движения спутника. Ось ординат y' дополняет систему до правой. В дальнейшем эту систему координат будем называть истинной звездной.
Для этой же цели широко применяется средняя звездная система координат (oxyz), начало которой совмещено с центром масс Земли. Ось аппликат z совпадает с ее средней осью вращения. Ось абсцисс x направлена в среднюю точку весеннего равноденствия, которая лежит в плоскости среднего экватора в эпоху t0. Ось ординат y дополняет систему до правой. В дальнейшем эту систему координат будем называть средней звездной.
В общем случае звездные системы координат не являются инерциальными, так как ось вращения Земли постоянно меняет свою ориентировку в пространстве.
Это сложное перемещение раскладывается на две составляющие: прецессию и нутацию. С прецессией связано пространственное перемещение так называемой средней оси вращения Земли (среднего полюса мира Pm), с нутацией – истинной (истинного полюса мира Pm'). Поэтому для соблюдения условия инерциальности задание пространственной ориентировки координатных осей на эпоху t0 является обязательным требованием.
В качестве некоторых аналогов инерциальной системы координат, пригодных для описания пространственного движения КА, Г. Вейсом предложены небесная и орбитальная системы координат [136]. В данной работе будет применяться прямоугольная геоцентрическая орбитальная система координат (О), начало которой совмещено с центром масс Земли. Ось аппликат совпадает с вектором кинетического момента. Ось абсцисс лежит в плоскости орбиты КА и образует с радиус-вектором r спутника угол, равный сумме трех углов: долготы восходящего узла, аргумента перигея и истинной аномалии. Ось ординат дополняет систему до правой [139]. Авторами работы [193] используется аналогичная система координат, называемая ими прямой равноденственной. Эта система координат является инерциальной.
В качестве другой орбитальной системы координат широко используется прямоугольная геоцентрическая орбитальная система координат (О'''), начало которой также совмещено с центром масс Земли. Ось аппликат ' совпадает с вектором кинетического момента. Ось абсцисс ' лежит в плоскости орбиты и направлена на спутник. Ось ординат ' находится в плоскости орбиты перпендикулярно к радиус-вектору КА и направлена в сторону его движения.
Данная система координат не будет инерциальной, так как оси ' и ' не сохраняют своей пространственной ориентировки.
Наряду с описанными выше системами в космической геодезии широко используются прямоугольные топоцентрические истинная (o x y z ) и средняя (o x y z ) системы координат, начала которых совпадают с точкой наблюдения o.
Координатные оси ориентируются в пространстве параллельно координатным осям соответствующих геоцентрических звездных систем.
определением звездного времени, фиксацией моментов наблюдений КА и т.п., используются различные системы времени. Одной из них является Всемирное определенного из обработки астрономических наблюдений суточных движений звезд, поправкой за движение истинного полюса Земли [1, 127 143] где tUT 1 и tUT 0 – моменты времени в шкалах UT1 и UT0 соответственно;
и – широта и долгота пункта наблюдения.
Для построения траекторий движения НП, искусственных и естественных небесных тел, для изучения неравномерности вращения Земли требуется некоторая шкала RT равномерного времени. В качестве таковой может быть принята шкала Атомного времени, не связанная с вращением Земли. Эта шкала основана на применении высокостабильных атомных эталонов частоты с относительной нестабильностью частоты порядка 3·10-16 (уход часов на 1 секунду происходит примерно за 70 миллионов лет) [27]. Такой цезиевый эталон времени создан в Национальном институте стандартов и технологии США. Поскольку Атомное время с высокой точностью воспроизводит шкалу RT, то его можно рассматривать как равномерное.
Существует несколько связанных между собой шкал Атомного времени, различающихся одна от другой выбором нуль-пункта [1, 135]:
ТА(SU) – шкала Атомного времени ГСВЧ (Россия);
TAS – ступенчатая шкала Атомного времени;
TAI – шкала Международного Атомного времени;
A.1– шкала Атомного времени, предложенная Морской обсерваторией США.
В качестве шкалы RT может также применяться шкала Всемирного координированного времени UTC, введенная с 1 января 1972 года для согласования систем Всемирного и Атомного времени. Поправка публикуется в бюллетене "Всемирное время и координаты полюса" серии Е. Здесь же публикуются разности (TAI-UTC), (TAI-TA(SU)) и (TA(SU)-UTC), позволяющие осуществить переход из одной шкалы времени в другую.
Поскольку вариации вращения Земли не подчиняются какому-либо известному в настоящее время закону, то и не существует аналитической зависимости между шкалой Всемирного времени UT1 и шкалой RT равномерного времени, необходимой для построения пространственных траекторий движения НП, естественных и искусственных небесных тел и т.д. Поэтому разность определяется либо после обработки результатов траекторных измерений КА, либо из астрономических наблюдений и также периодически публикуется в бюллетенях.
1.3 Преобразования координат В космической геодезии невозможно обойтись без разнообразных координатных преобразований [76]. Это обусловлено тем, что пространственные положения различных объектов и рассматриваемых явлений определены в различных системах координат, а задача оценивания должна выполняться в единой и удобной для работы системе координат. Ниже рассмотрим основные координатные преобразования, основой для которых служат матрицы простых вращений. Они применяются для преобразования векторов из одной прямоугольной системы координат в другую. Это матрицы: R1() – при вращении вокруг оси “ox”, R2() – при вращении вокруг оси “oy”, R3() – при вращении вокруг оси “oz”. Положительным углом вращения считается такой угол, при котором вращение системы координат происходит против часовой стрелки, если смотреть с положительного конца оси вращения.
Преобразование средних земных координат (OXYZ)G в истинные земные осуществляется двумя последовательными поворотами системы (OX'Y'Z')G (OXYZ)G: сначала относительно оси OXG на угол yP, затем относительно оси OYG (уже совпадающей с осью OY) на угол xP (рисунок 1.1). Здесь величины xP и yP являются координатами истинного (мгновенного) полюса Р' Земли относительно среднего полюса Р.
PG – средний меридиан Гринвича, P'G' – истинный меридиан Гринвича, QQ – средний экватор, Q'Q' - истинный экватор, 0P – средняя ось вращения Земли, 0P' – истинная ось вращения Земли, G – точка пересечения среднего меридиана Гринвича со средним экватором, G' – точка пересечения истинного меридиана Рисунок 1.1 – Связь истинных и средних земных систем координат Матричное преобразование средних земных координат в истинные может быть представлено в виде:
а истинных земных координат в средние – где RG [XYZ ]T – геоцентрический радиус-вектор НП в общеземной системе координат;
RG [X 'Y 'Z ' ]T – геоцентрический радиус-вектор НП в истинной земной системе координат.
Пусть положение НП определено в референцной системе координат (O''XYZ)Г радиус-вектором R [XYZ ]T, а в общеземной (OXYZ)G – радиус-вектором RG [XYZ ]T (рисунок 1.2).
(OXYZ)'' – промежуточная (квазигеоцентрическая) система с осями координат параллельными Рисунок 1.2 – Связь геоцентрических и квазигеоцентрических выполняется с помощью матричного выражения [8, 22, 89, 95] где dR [dX dY dZ ]T – трехмерный вектор-столбец смещения начала О'' системы координат (O''XYZ)Г относительно начала О системы (OXYZ)G;
[ X Y Z ]T – трехмерный вектор-столбец малых углов поворота координатных осей системы координат (O''XYZ)Г относительно осей системы (OXYZ)G;
(O''XYZ)Г.
Обратное преобразование имеет вид:
В формулах (1.51) и (1.52) матрица малых поворотов координатных осей перемножения трех матриц вращения, имеет вид [136]:
Компоненты векторов dR, и поправка k называются элементами взаимного ориентирования систем координат.
После выполнения простых матричных преобразований, выражения (1.51) и (1.52) могут быть записаны в следующем виде:
Здесь матрицы D1 и D2 имеют вид Если пространственное положение НП задано в общеземной системе координат (0XYZ)G радиус-вектором RG [XYZ ]T, а в горизонтной (0XYZ)H задано радиус-вектором R H [XYZ ]T (рисунок 1.3), то преобразование этого вектора из общеземной системы в горизонтную осуществляется сначала поворотом системы (OXYZ)G относительно оси OZG на угол (90+L), а затем относительно оси OXG на угол (90- B):
Рисунок 1.3 – Связь горизонтной и общеземной систем координат После перемножения матриц поворота преобразование (1.57) будет иметь вид:
В нем Обратное преобразование имеет вид Элементы матрицы Q являются функциями прямоугольных земных координат НП и рассчитываются по формулам:
sinB = {1+e2·(1-ZG2/RG2)}· ZG/ RG, cosB = 1-sin 2 B, где B, L – геодезические широта и долгота НП соответственно;
e – эксцентриситет земного эллипсоида.
Преобразование вектора RG положения НП из общеземной системы координат в звездную описывается выражением где R и R – геоцентрические радиус-вектор положения и вектор скорости НП в звездной системе координат;
P и P – матрицы преобразования.
Выражение (1.61) является кинематической моделью движения НП в пространстве.
Вид матриц P и P в выражении (1.61) определяется выбранной координатной системой, которая на интервале времени T анализа состояния нелинейной ДС будет удовлетворять требованиям инерциальной системы. Матричное преобразование (1.61) должно быть простым и, одновременно, обеспечивать точность, соответствующую современному уровню качества траекторных измерений КА.
ортогональных матриц вращения, учитывающих вращение Земли вокруг своей оси, колебания этой оси вращения, прецессию по прямому восхождению, и склонению на отрезке времени t = t - t0 (t0 – начальный момент времени, t – текущий момент времени, соответствующий очередному измерению), а также нутацию по прямому восхождению t и t 0, склонению t и t 0, наклону t и t 0 в эпохи t и t0 соответственно [140]:
Здесь обозначено рассчитывается по формуле [140] Звездное время S0 определяется как где JD(d 0 ) – юлианский день 0 UT1 даты d0 задания НУ движения КА;
JD(d ) – юлианский день 0 UT1 текущей даты d;
– общая прецессия по прямому восхождению истинной точки весеннего равноденствия за период t0 - tБ (tБ – начало Бесселева тропического года 1950.0) [135].
В формулах (1.63) – (1.67) моменты времени t и t0 должны быть заданы в шкале Всемирного времени UT1.
Величины,, x p, y p и UT могут быть заранее аппроксимированы на отрезке времени [t0, t] степенным полиномом с аргументом t = t - t0, который задается в равномерной шкале времени RT [136, 140]:
где k, k, xk, yk,UTk – коэффициенты полиномов;
Матрицу P получаем дифференцированием матрицы P по времени [49, 53]:
где – скорость изменения нутации по наклону на отрезке времени t.
В выражении (1.71) величина определяется по формуле где и – скорости изменений прецессии и нутации по склонению на отрезке времени t.
Выражения для вычисления и получаем дифференцированием по времени полиномов, определяющих значения параметров и :
Выражения (1.69) и (1.70) представляют собой кинематические модели изменения ПВЗ, а величины x0,.., xL1, y0,..., y L1, UT0,..., UTL1 являются параметрами этих моделей.
Для расчета прямоугольных координат векторов положения и скорости спутника по его регулярным элементам (см. 1.5) удобнее использовать (рисунок 1.4). Ее образует ортонормированный базис начало которого совмещено с центром масс Земли. Компоненты базисных векторов в выражении (1.73) являются функциями регулярных элементов орбиты спутника и определяются по формулам (В.4) (приложение В).
i – наклонение орбиты, – долгота восходящего узла, – аргумент перигея, v – истинная аномалия, r – геоцентрический радиус-вектор КА, р – перигей орбиты, (oxyz) – звездная Рисунок 1.4 – Инерциальная орбитальная система координат (O) В этой системе координат ось абсцисс лежит в плоскости орбиты и образует с радиус-вектором r спутника угол l = + + v. Согласно правилам линейной алгебры матричное преобразование радиуса-вектора положения r и скорости r спутника из звездной системы координат (oxyz) в орбитальную (O) будет иметь вид [101]:
где re=[ ]T – радиус-вектор спутника в орбитальной системе координат;
r = [x y z]T – радиус-вектор спутника в звездной системе координат;
re = [ ]T – вектор геоцентрической скорости спутника в орбитальной системе координат;
r [x y z ]T – вектор геоцентрической скорости спутника в звездной системе координат.
Выполнив скалярное произведение векторов в выражении (1.74), получим проекции вектора r на оси орбитальной системы координат [171]:
Дифференцирование по времени правой и левой частей равенства (1.76) позволяет получить зависимости для определения проекций вектора скорости КА на оси орбитальной системы координат Для практического использования в плане расчета величин и следует расстояния r и скорости r спутника в функциях от регулярных элементов орбиты определяются по формулам (1.138) и (В.5) (приложение В) соответственно.
Обратные к выражениям (1.74) и (1.75) матричные преобразования имеют вид или Наряду с описанной выше системой, в космической геодезии также активно используется и орбитальная система координат (O''') (рисунок 1.5).
u – аргумент широты КА, (oxyz) – звездная система координат Рисунок 1.5 – Орбитальная система координат (O''') Она удобна для пересчета параметров движения КА (например, кеплеровых элементов) в прямоугольные координаты векторов положения и скорости спутника, для учета возмущающих ускорений, влияющих на движение спутника по орбите и т.д. У нее ось абсцисс O' постоянно направлена на КА, ось O' совпадает с вектором кинетического момента, ось O' лежит в плоскости орбиты и ортонормированный базис S, T и W, задающий направление координатных осей ', ' и '. Начало этой системы совмещено с центром масс Земли.
Матричные преобразования векторов положения r и скорости r спутника из орбитальной системы координат (O''') в звездную (oxyz) имеют вид [144] где r – геоцентрическое расстояние до КА;
Vr – радиальная скорость КА;
Vn – трансверсальная скорость КА.
В выражениях (1.82) и (1.83) скорости Vr и Vn являются функциями кеплеровых элементов орбиты. Они вычисляются согласно выражениям [8, 144]:
где – гравитационный параметр Земли;
e – эксцентриситет орбиты КА;
v – истинная аномалия;
p – фокальный параметр.
Базисные векторы S, T и W также являются функциями кеплеровых элементов. Они определяются по формулам [89, 144, 145]:
где i – наклонение орбиты КА;
– долгота восходящего узла орбиты;
u – аргумент широты КА.
1.4 Общие положения и математическая постановка задачи определения геодезических параметров по результатам радиотехнических траекторных измерений спутников Для оценивания состояния нелинейной динамической системы, состоящей из нескольких КА и сети НП, в космической геодезии используются несколько методов, в основном различающихся между собой набором компонент вектора оцениваемых параметров и применяемой теории построения спутниковых орбит.
К таким методам относятся: динамический, полудинамический [8, 136, 144, 178], коротких дуг [19, 146, 166, 172], транслокации [112, 172, 182, 190], точечного позиционирования [166, 172, 182, 190]. К основным задачам, которые позволяют решать эти методы, относятся:
– уточнение параметров гравитационного поля Земли;
– определение пространственных положений НП в заданной геодезической системе координат и создание на поверхности Земли единой мировой опорной геодезической сети;
– определение ЭВО геодезических систем координат;
– определение ПВЗ с целью изучения движения земных полюсов и неравномерности вращения 3емли вокруг своей оси.
динамической системы сводится к формированию и решению уравнений измерений, использующих различные виды траекторных измерений КА. К ним относятся: оптические (фотографические, лазерные), запросные и беззапросные радиотехнические (доплеровские, радиодальномерные). Все эти виды измерений различаются между собой по точности, оперативности, универсальности (в смысле их всепогодного и круглосуточного применения).
Применительно к орбитальному методу космической геодезии далее рассмотрим теорию формирования уравнений поправок для анализа поведения принятой нелинейной ДС на заданном отрезке времени T по результатам беззапросных радиотехнических траекторных измерений КА [86, 133]. Заметим, что эти уравнения являются составной частью уравнений поправок динамического метода. В качестве вектора q параметров состояния этой системы в текущий момент времени t примем вектор-столбец где Yi – шестимерный вектор-столбец элементов орбиты i - го КА (i = 1, 2,..., g;
g – число КА);
Rj, R j – определяемые из выражения (1.61) трехмерные векторы-столбцы прямоугольных координат положения и скорости j- го НП в звездной системе координат (j = 1, 2,..., l; l – число НП, результаты измерений с которых приняты в обработку);
dR – трехмерный вектор-столбец смещения начала референцной системы координат относительно начала общеземной системы;
– трехмерный вектор-столбец малых углов поворота координатных осей референцной системы координат относительно осей общеземной системы;
k – поправка к масштабу референцной системы координат;
– трехмерный вектор-столбец ПВЗ вида = [x p y p UT ]T.
Будем считать, что в результатах измерений присутствуют сингулярные ошибки вида (1.26) и (1.38). Поэтому задачу оценивания сведем к отысканию компонент вектора q0 начального состояния системы в некоторый заданный момент времени t0 и компонент вектора мешающих параметров.
Уравнения измерений (1.5) являются нелинейными функциями векторов q и окрестности априорных значений векторов q и Г, удерживая при этом только линейные члены разложения. В результате можно записать где q, Г – априорные значения векторов q и Г соответственно.
В уравнениях (1.87) разности q и Г определяются как Рассматриваемая динамическая система является прогнозируемой, т.е. для нее выполняется равенство (1.2). Следовательно, справедливо дифференциальное соотношение в котором матрица частных производных рассчитывается для априорных значений вектора q'0 НУ, а вектор q0 находится как С учетом выражений (1.20) и (1.89), равенство (1.87) примет вид:
Здесь матрица П рассчитывается по формуле (1.20) с использованием априорно заданного вектора q' параметров состояния анализируемой ДС.
Подставив равенство (1.91) в выражение (1.15), с учетом зависимости (1.5), в общем виде получим систему уравнений поправок на текущий момент времени t где q0 – r-мерный вектор-столбец поправок к вектору q'0 априорных значений вектора q0 параметров состояния ДС на начальный момент времени t0;
Z – d-мерный вектор-столбец результатов измерений вида (1.8);
Z – d-мерный вектор-столбец первичных измеряемых псевдопараметров использованием векторов q'0 и Г';
Vz – d-мерный вектор-столбец поправок к вектору Z.
В уравнениях (1.92) векторы Z и Vz определяются по формулам где – случайный d-мерный вектор-столбец ошибок измерений;
– случайный d-мерный вектор-столбец ошибок моделей измерений.
Так как рассматриваемая ДС на отрезке времени Т является прогнозируемой, то задача оценивания сводится к отысканию компонент расширенного вектора Q начального состояния системы на некоторый начальный момент времени t0.
Вектор Q0 сформирован на основе вектора q из равенства (1.86), в котором согласно матричному выражению (1.61) вместо векторов R и R принят вектор RG положения НП в общеземной системе координат. В результате расширенный вектор оцениваемых параметров для одного сеанса измерений примет вид где Y0i – шестимерный вектор-столбец НУ движения i-го КА (i = 1, 2,..., g;
g - число КА);
RGj – трехмерный вектор-столбец прямоугольных координат j-го НП (j = 1, 2,..., l; l – число определяемых НП) в общеземной системе;
Г – вектор-столбец мешающих параметров в одном сеансе измерений;
0 – вектор-столбец параметров кинематических моделей изменения ПВЗ.
Вектор-столбец 0 в равенстве (1.93) состоит из коэффициентов полиномов, представленных выражениями (1.69) и (1.70), и имеет вид Теперь система уравнений поправок (1.92) для совместного оценивания НУ движения КА, ПВЗ и координат НП, заданных в общеземной системе, примет вид где Г, Y0, RG и 0 – векторы-столбцы поправок к одноименным априорным (приближенным) значениям векторов-столбцов, заданных на начальный момент времени t0.
Если спутниковые измерения выполнены с НП, пространственные положения которых известны с достаточно высокой точностью в референцной системе координат (OXYZ)Г, то присутствие вектора RG в уравнениях (1.95) будет обусловлено несовпадением референцной и общеземной систем координат. Этот вектор легко получается из формулы (1.54) как разность Подставив значение вектора RG из равенства (1.96) в выражение (1.95), получим общий вид системы уравнений поправок для совместного оценивания по результатам измерений одного сеанса начальных условий движения КА, ПВЗ и ЭВО референцной системы координат относительно общеземной:
Теперь получим общий вид системы уравнений поправок для вторичных динамической системы по результатам измерений одного сеанса. Для этого воспользуемся зависимостью (1.29), в которой вектор-столбец вторичных измеряемых псевдопараметров представлен формулами (1.30) и (1.31).
Очевидно, что тогда системы уравнений поправок, аналогичных по назначению уравнениям (1.95) и (1.97), имеют вид [33, 34]:
где E – матрица частных производных вида (1.32), являющаяся матрицей коэффициентов перед неизвестным вектором Г, рассчитывается по априорным (приближенным) значениям вектора q = q.
Перемножив левые и правые части уравнений (1.98) и (1.99) на матрицу, получим общий вид системы уравнений поправок для вторичных измеряемых псевдопараметров при совместном оценивании начальных условий движения КА, ПВЗ и координат НП а также систему уравнений поправок для совместного оценивания НУ движения КА, ПВЗ и ЭВО где – d-мерный вектор-столбец результатов измерений вторичных измеряемых параметров;
– d-мерный вектор-столбец вторичных измеряемых псевдопараметров;
V – d-мерный вектор-столбец поправок к вектору.
В уравнениях (1.100) и (1.101) векторы-столбцы и определяются по формулам Ниже рассмотрим широко используемый метод статистической обработки результатов траекторных измерений КА, в котором за критерий качества принят статистический обработки измерений называется методом наименьших квадратов.
Пусть в течение нескольких сеансов на отрезке времени T в дискретные моменты t1, t2, …, tm (m – общее число измерений, взятых в обработку) выполнены неравноточные траекторные измерения КА, которым соответствует модель измерений (1.5). Будем считать, что вектор q = [q(t1) q(t2) … q(tm)]T параметров состояния принятой нелинейной ДС на эти моменты времени определяется по известным начальным условиям q 0 равенства (1.3) с помощью моделей движения для КА – формулы (1.137), для НП – формулы (1.61), для ПВЗ – формулы (1.69) и (1.70). В этом случае расширенный n-мерный вектор оцениваемых параметров, определяемый выражением (1.93), примет вид Здесь совокупный вектор-столбец мешающих параметров размерности k·p и размерность n вектора Q0 определяются как где k – размерность вектора Г;
p – число сеансов измерений;
g – число КА;
l – число определяемых НП;
L1 – порядок полиномов (1.69) и (1.70).
Будем также считать, что математическое ожидание совместного влияния случайных ошибок измерений =[(t1) (t2) … (tm)]T и модели измерений =[(t1) (t2) … (tm)]T равно нулю M( + ) = 0, а ковариационная матрица оцениваемых параметров Q0, для которого квадратичная форма достигает своего минимума. Здесь вектор-столбец вторичных измеряемых псевдопараметров U(Q0 ), рассчитываемый по приближенным значениям вектора Q 0 оцениваемых параметров и результатам измерений вторичных измеряемых параметров, имеет вид где р – общее число сеансов;
m – общее число измерений.
Для выполнения этого требования должно соблюдаться равенство где q, – оценки векторов q и Г.
Согласно правилу дифференцирования векторных функций получим Последнее выражение представляет собой систему нелинейных уравнений. В ней функциональная матрица есть прямоугольная mxn матрица частных производных от измеряемых векторфункций по компонентам вектора Q 0.
Система нелинейных уравнений (1.112) решается итерационными методами.
Наибольшее распространение получил метод последовательных приближений Гаусса-Ньютона [87, 88, 144, 155]. Использование этого метода предусматривает линеаризацию измеряемых вектор-функций в окрестности априорных (приближенных) значений оцениваемых параметров. Применив метод Гаусса-Ньютона для решения системы нелинейных уравнений (1.112), получим выражение В нем вектор-столбец поправок к приближенным значениям вектора оцениваемых параметров Q (0 j ) равенства (1.105) определяется как Выражение (1.114) позволяет найти уточненную оценку Q(j+1) вектора оцениваемых Итерационный процесс заканчивается при выполнении критерия где j – номер приближения;
– вектор-столбец априорно заданных малых параметров, размерность которого соответствует размерности вектора оцениваемых параметров равенства (1.105).
полученных оценок Q 0 вычисляется по формуле Используя общепринятые обозначения, систему уравнений поправок вида (1.101) для всех видов измерений, принятых в обработку, можно записать как где X – вектор-столбец неизвестных вида (1.115), состоящий из поправок к приближенному значению вектора оцениваемых параметров Q 0 равенства (1.105);
V – вектор-столбец поправок к вектору результатов измерений;
F – m-мерный вектор-столбец правой части.
Система уравнений (1.119) представляет собой несовместную СЛАУ. В ней вектор-столбец F определяется по формуле рассчитываемый по приближенным значениям вектора оцениваемых параметров, имеет вид Как уже отмечалось выше, решение нелинейной задачи МНК включает процедуру замены нелинейной задачи оценивания на линейную, т.е. переход от нелинейных моделей измерений к линейным. В частности, решение системы нелинейных уравнений методом Гаусса-Ньютона предусматривает последовательное выполнение ряда линейных задач до достижения принятого критерия окончания итеративного процесса. Иными словами, неоднократно решается линейная задача наименьших квадратов, суть которой заключается в следующем [15, 102, 148]. Задана несовместная СЛАУ(1.119). Требуется найти такой действительный вектор X, для которого евклидова норма • Е вектора A X F будет минимальной:
где X – оценка вектора X (или вектор решений).
По отношению к несовместной системе уравнений V будет вектором невязок [5], значение которого определяется равенством Геометрически условие (1.122) означает, что необходимо найти такую точку с координатами A X в n системе m-мерных базисных векторов которая ближе всего находится к вектору F по евклидовой норме. Этому требованию удовлетворяет точка, для которой вектор невязок V будет ортогонален каждому базисному вектору ai (i = 1,..., n), то есть Матричная запись (1.126) представляет собой совместную СЛАУ известную под названием линейной системы нормальных уравнений с квадратной nxn матрицей коэффициентов AT·A и вектором правой части AT·F. Решением этой системы является вектор Обобщенным решением (псевдорешением) системы (1.123), в смысле наименьших квадратов, будет любой вектор X, для которого евклидова норма вектора невязок V из (1.124) достигает наименьшего значения (1.122). Решение X, имеющее наименьшую евклидову норму min X, будет нормальным обобщенным решением (нормальным псевдорешением) [115]. Обобщенные решения и только они удовлетворяют СЛАУ (1.126). Если матрица A неполного ранга, то существует бесчисленное множество обобщенных решений. Среди них имеется нормальное псевдорешение, которое всегда существует и единственно.
Как будет показано ниже, единственное решение СЛАУ (1.123) находится из выражения [87, 102, 115] где A – псевдообратная (либо эффективная псевдообратная) матрица для матрицы A.
На заключительном этапе линейной задачи МНК можно находить вектор X из решения СЛАУ с квадратной (равенство (1.127)) либо прямоугольной (равенство (1.128)) матрицей коэффициентов. Зачастую такие матрицы определены с ошибками. Кроме того, элементы этих матриц в оперативной памяти компьютера представляются также с ошибками. Вполне возможно, что матрица коэффициентов будет неполного ранга (факт нерегулярной постановки измерительной задачи) или близка к этому (факт некорректной постановки измерительной задачи). А это приводит либо к неопределенности решения, либо к повышению чувствительности решения СЛАУ к ошибкам входных данных. В связи с перечисленными выше положениями, на этапе постановки задачи оценивания необходимо выбирать оптимальную стратегию ее решения и правильно интерпретировать полученные результаты. Поэтому необходимо применять современные методы вычислительной алгебры, чтобы обеспечить регулярную и корректную постановку задачи оценивания расширенного вектора параметров состояния конкретной нелинейной ДС [51, 80].
Заключительный этап выполнения линейной задачи МНК – решение СЛАУ (1.119) представим в виде решающей подсистемы и рассмотрим методику исследования ее устойчивости. У этой системы входом служат матрица коэффициентов и вектор правой части, а выходом – векторы X и V.
По определению А. Н. Тихонова и В. Я. Арсенина решение СЛАУ будет устойчивым, если малым изменениям входных данных соответствуют малые изменения вектора решения [141]. Решение неустойчиво, если малые изменения входных данных могут привести к большим изменениям вектора решения.
Линейные системы, чувствительные к изменению входных данных, называются плохо обусловленными или неустойчивыми [23, 115]. Такие СЛАУ содержат плохо обусловленные матрицы коэффициентов [152].
При наличии возмущений A и F матрицы коэффициентов и вектора правой части фактически решается система уравнений где X – n -мерный вектор-столбец решения этой возмущенной системы.
Верхний предел его относительной ошибки согласно [22, 102] определяется неравенством с параметрами
A E X AX A E X AX E
где (А) – число обусловленности матрицы A.Из неравенства (1.131) следует, что устойчивость решения СЛАУ к изменению входных данных существенно зависит от числа обусловленности (А). Оно выступает в роли коэффициента увеличения совместного влияния относительных ошибок матрицы коэффициентов А и вектора правой части F на результат: с увеличением числа обусловленности, соответственно, увеличивается и относительная ошибка решения X. Если число обусловленности мало (близко к единице), то матрица A будет хорошо обусловленной. Если число обусловленности велико, то данная матрица будет плохо обусловленной [115, 148, 152].
Из анализа выражений (1.131) и (1.132) следует, что на число итераций существенное влияние оказывают ошибки задания матрицы коэффициентов A и вектора F. Ошибка расчета вектора F зависит от точности взятых в обработку измерений и точности математических моделей этих измерений. Если матрица A плохо обусловлена ((А) 1), то необходимо принять меры к повышению точности расчета этой матрицы и вектора правой части так, чтобы выполнялось условие (А)·А 1. В этом случае соблюдается неравенство X 1 и итеративный процесс решения системы нелинейных уравнений будет сходящимся. Кроме того, зная ошибку А до решения СЛАУ и применяя процедуру сингулярного разложения (см. приложение Б), можно определить эффективный ранг матрицы коэффициентов. Последнее позволяет уверенно найти нормальное обобщенное решение заданной системы уравнений.
Переходя к методам решения линейных систем, отметим, что выражение (1.128) здесь и в целом ряде работ используется только для краткости записи.
Обычно при решении СЛАУ псевдообратная (либо обратная) матрица не вычисляется, так как более эффективно находить решение с применением методов разложения матрицы коэффициентов на матрицы-сомножители (обычно на две или три). Для квадратных матриц чаще всего используются разложения Холесского, Грама-Шмидта, Гаусса, а для прямоугольных – Хаусхолдера, SVD (сингулярное разложение, см. приложение Б). Достоинство этих методов разложения состоит в том, что они позволяют путем ортогональных преобразований привести исходную матрицу к более простому виду: треугольному, диагональному и т.д. При этом ортогональные преобразования не вносят дополнительных ошибок в вычисления.
Упомянутые методы разложений достаточно хорошо описаны в специальной литературе, посвященной вычислительным аспектам линейной алгебры.
В заключение отметим, что при решении системы линейных уравнений поправок лучше применять метод SVD по двум основным причинам (см. приложение Б). Во-первых, здесь не требуется приводить матрицу коэффициентов к квадратному виду, что не увеличивает ее обусловленности. Во-вторых, этот метод позволяет определять количественные условия наблюдаемости и, не привлекая дополнительных средств анализа, учитывать их на этапе решения уравнений.
1.5 Выбор системы элементов орбиты спутника и фундаментальной матрицы Другим актуальным вопросом для обеспечения регулярной постановки задачи оценивания является выбор оцениваемых параметров Y орбиты спутника и формирование на их основе фундаментальной матрицы – составной части матрицы коэффициентов уравнений поправок.
Если в качестве компонент вектора Y выбраны кеплеровы элементы, то для орбит с малыми эксцентриситетами (e < 0.003) матрица коэффициентов системы уравнений поправок становится плохо обусловленной [8]. Кроме того, дифференциальные уравнения движения вида (1.3) в кеплеровых элементах имеют особенности для круговых, полярных и экваториальных орбит. По этим причинам целесообразнее применять прямоугольные координаты векторов положения r = [x y z]T и скорости r =[x y z ]T КА, канонические элементы орбиты, а также различные модификации кеплеровых элементов [8, 122, 144].
В ряде работ [171, 173, 193] предложена "несингулярная" система элементов орбиты, являющаяся функцией кеплеровых и не имеющая особенностей для круговых и полярных орбит:
где a, e – большая полуось и эксцентриситет орбиты;
i – наклонение орбиты;
– долгота восходящего узла;
– аргумент перигея орбиты;
M – средняя аномалия;
– средняя долгота спутника.
Авторами работ [91, 138] была предложена аналогичная несингулярной системе (1.134) система регулярных элементов. Отличие заключается только в том, что у регулярных элементов вместо большой полуоси и средней аномалии используется фокальный параметр р и истинная аномалия v. Это позволило получить выигрыш по времени при интегрировании дифференциальных уравнений движения. По этой причине в данной работе за основную принята система регулярных элементов орбиты, образующих шестимерный векторстолбец Компоненты этого вектора рассчитываются по формулам оскулирующих регулярных элементов (1.136) не имеет особенностей в случае круговых, полярных и экваториальных орбит [137, 140]. Эта система имеет вид:
В них где – гравитационный параметр Земли;
r – геоцентрическое расстояние до КА [171];
S, Т, W – компоненты вектора возмущающего ускорения в подвижной орбитальной системе координат (S,T,W) (см. рисунок 1.5).
Здесь и далее точка над переменной означает дифференцирование параметра по времени. Система дифференциальных уравнений движения для несингулярных элементов орбиты (1.134) приведена в работе [173].
Одной из трудоемких операций при формировании матрицы коэффициентов в уравнениях поправок является вычисление матрицы частных производных от измеряемых вектор-функций по начальным условиям, являющейся матрицей коэффициентов перед вектором поправок Y0 к регулярным элементам орбиты (1.135) в уравнениях (1.95), (1.97), (1.100) и (1.101). Она может быть представлена в виде произведения двух матриц:
Здесь есть шестимерный вектор-столбец, состоящий из прямоугольных координат векторов положения r и скорости r КА в текущий момент времени t T. Первая матрица-сомножитель (градиентная матрица первого рода [13, 88, 159]) в равенстве (1.139) является прямоугольной mх6 матрицей частных производных от измеряемых вектор-функций зависимости (1.30) по текущим значениям компонент векторов r и r КА Вторая матрица-сомножитель (фундаментальная матрица первого рода [88, 159]) в равенстве (1.139) является квадратной 6х6 матрицей частных производных от текущих значений компонент векторов r и r (1.140) по начальным условиям движения КА в регулярных элементах Градиентная матрица отражает специфические свойства каждого вида измерений. В частности, она показывает, как изменяются измеряемые параметры с изменением пространственного положения и скорости КА. Фундаментальная матрица является общей для всех видов траекторных измерений и определяет общие свойства решения системы линейных уравнений поправок.
фундаментальной матрицы является сложной самостоятельной задачей. Вид этой матрицы изменяется в зависимости от выбранной системы элементов орбиты. В работах [8, 98, 122] приводятся аналитические выражения для расчета фундаментальной матрицы в кеплеровых элементах, а в работе [15] – в прямоугольных координатах векторов r и r.
представить в виде произведения двух матриц Первым сомножителем в выражении (1.143) является квадратная 6х6 матрица частных производных от текущих значений прямоугольных координат X вида (1.140) векторов положения и скорости КА по текущим значениям регулярных элементов Y его орбиты Вторым сомножителем является квадратная 6x6 матрица частных производных от текущих значений регулярных элементов орбиты по их начальным значениям (матрица изохронных производных) Вывод аналитических выражений для расчета элементов этих матриц представлен в приложении В и в приложении Г. Анализ точности вычисления матрицы изохронных производных приведен в приложении Д.
Общий вид фундаментальной матрицы Ф можно получить, перемножив матрицы в равенстве (1.143). При этом следует учитывать зависимости (В.1) (приложение В) для матрицы С и (Г.14) (приложение Г) для матрицы изохронных производных С0. В связи с этим можно записать или Для нахождения аналитических выражений элементов матрицы Ф воспользуемся формулами (В.12), (В.13) (приложение В) и (Г.41) (приложение Г) расчета элементов матриц C и С0. Фундаментальную матрицу получим с допущением p = p0.
Для элементов первого и шестого столбцов можно записать где t0 – начальный момент времени;
t – текущий момент времени.
Для элементов четвертого и пятого столбцов фундаментальной матрицы справедливы следующие равенства:
Далее получим аналитические выражения для расчета элементов второго столбца фундаментальной матрицы. Можно записать, что эксцентриситета e орбиты то с учетом зависимостей (1.136) и (Г.40) для величин g, h и имеем:
Откуда запишем или или Получим аналитические выражения для расчета элементов третьего столбца фундаментальной матрицы. Можно записать или или Полученные производные для расчета элементов фундаментальной матрицы Ф согласно [36, 44] можно представить в виде:
Параметры K, L, H, и находятся из выражений где – гравитационный параметр Земли;
r – геоцентрическое расстояние до КА в текущий момент времени t;
r0 – геоцентрическое расстояние до КА в начальный момент времени t0.
1.6 Оценка влияния изменения гравитационного поля Земли на движение космического аппарата сопротивления атмосферы, лунно-солнечного притяжения, прямого светового давления от Солнца и диффузного отражения света от Земли, оказывают свое влияние целый ряд других возмущений. Некоторые из них рассмотрены, например, в работе [2]. Кроме этого, в работе [129] говорится о том, что влияние крупных водохранилищ и изменение масс озер могут привести к изменению в коэффициентах низких гармоник модели геопотенциала уже в третьей значащей цифре. В связи с этим в данном подразделе рассматривается аналитическое решение задачи оценки влияния изменения параметров модели ГПЗ на движение КА в пределах неуправляемого участка его полета. Работа сделана с участием соискателя в рамках выполнения фундаментальной НИР “Исследование динамической составляющей гравитационного потенциала Земли и его характеристик, построение их математических моделей” [45, 92], которая стала продолжением исследований в рамках фундаментальной НИР “Разработка методов изучения динамики гравитационного поля и фигуры Земли планетарного, регионального и локального характера” [16 – 18, 130, 131].
С позиций системного подхода задачу оценивания изменения ГПЗ на движение КА можно разделить на две подзадачи: первая – определение изменения ГПЗ, обусловленного геодинамическими, техногенными и другими процессами, а также ошибочностью определения этого поля [45].; вторая – оценка влияния этих изменений на движение КА [83]. Первая подзадача сводится к определению изменения числовых значений коэффициентов сферических функций. Вторая подзадача – к определению влияния изменения этих коэффициентов на движение КА в околоземном пространстве. Ниже представим математический алгоритм решения второй подзадачи.
Для ее решения воспользуемся системой дифференциальных уравнений Э = [a e i M]T орбиты КА где S, T, W – нормированные значения возмущающего ускорения, разложенного в базисе S, T и W орбитальной системы координат (O''') (см. 1.3).
а – большая полуось орбиты спутника;
р – фокальный параметр орбиты;
формуле [13]:
где S, T, W – значения возмущающего ускорения, разложенного в ортонормированном базисе S, T и W орбитальной системы координат (O''');
Обозначив нормированные значения возмущающего ускорения как векторстолбец а значения возмущающего ускорения - как вектор-столбец равенство (1.162) можно представить в виде:
Правые части уравнений (1.161), являющиеся производными кеплеровых Будем считать, что на движение КА оказывает свое влияние только ГПЗ. Тогда, разложив правые части этих уравнений в ряд и ограничившись при этом только линейными членами разложения, можно записать:
где C, S – коэффициенты разложения ГПЗ в ряд по сферическим функциям;
n и m – порядок и степень разложения ГПЗ в ряд.
Первое слагаемое A0 в разложении (1.166) – это вектор-столбец правой части системы уравнений Ньютона-Лагранжа, полученных в стационарном ГПЗ, второе слагаемое В – вектор-столбец изменений правых частей уравнений движения (1.161), обусловленных изменением ГПЗ (или ошибками его представления).
Под вектором возмущающего ускорения R будем понимать градиент от возмущающего гравитационного потенциала V Земли где ms – масса спутника.
определяется как [2, 21, 89, 135] где Pnm(sin) – функции Лежандра;
RЭ – средний экваториальный радиус Земли;
Разложение вектора возмущающего ускорения R в средней звездной системе координат (oxyz) можно представить следующим образом:
где x, y, z – прямоугольные координаты КА в средней звездной системеж r – геоцентрическое расстояние до КА;
, – геоцентрические широта и долгота КА;
xG, yG, zG – прямоугольные координаты КА в общеземной системе.
Введя обозначения выражение (1.169) можно записать как Чтобы получить разложение вектора R в орбитальной системе координат (O'''), необходимо выполнить следующее матричное преобразование представить в виде произведения частных производных Для удобства дальнейших аналитических выкладок, систему уравнений (1.161) представим в матричном виде Здесь матрица D является функцией кеплеровых элементов орбиты КА. Она имеет вид воспользовавшись равенствами (1.165), (1.167), (1.172) и (1.174):
где Е – единичная 3x3 матрица.
Обозначив третью производную в матричном преобразовании (1.174) как Q и приняв во внимание выражения (1.85) для расчета компонент базисных векторов S, T и W, запишем Матрица Q является матрицей преобразования векторов из инерциальной системы координат (oxyz) в орбитальную (O''') [145].
Выражение для определения пятой производной в равенстве (1.173) найдем, воспользовавшись зависимостью (1.171):
Преобразование геоцентрического радиус-вектора r КА, заданного в звездной системе координат (oxyz), в геоцентрический радиус-вектор rG, заданного в земной системе координат (oxyz)G, представляет следующий вид:
где R3(S) – матрица вращения на угол S, являющимся звездным временем.
В этом случае матрица Р размерности 3x3 в равенстве (1.170) примет вид Для определения элементов матрицы О размерности 3x3 в равенстве (1.170) воспользуемся следующими формулами:
Выполнив их дифференцирование, получим согласно [145]:
Аналитические выражения для расчета элементов матрицы N размерности 1x в равенстве (1.170) можно найти, например, в работе [21]. Они выглядят следующим образом:
В зависимостях (1.184) частные производные от функций Лежандра Pnm(sin) по геоцентрической широте спутника определяются следующим образом [21]:
Ее элементы получим, последовательно дифференцируя матрицу (1.184) по коэффициентам C20, C21, C22, C30, C31, C32, C33, …, S20, S21, S22, S30, S31, S32, S33, …:
«