Содержание к диссертации
Введение 5
1 Краевые задачи теории потенциала и методы их решения 16
Потенциалы тяготения и силы тяжести 16
Типы краевых задач 20
Методы решения краевых задач 21
Общие сведения о приближенных методах решения 21
Метод взвешенных невязок 22
Метод коллокаций 24
Метод конечных разностей 25
Метод коллокаций с подобластями 28
1.3.6. Метод Бубнова-Галеркина 29
Метод наименьших квадратов 29
Слабая формулировка краевой задачи 30
Метод конечных элементов 30
Метод граничных интегральных уравнений 32
Решение уравнения Лапласа методом разделения переменных 3 8
Метод точечных масс 41
1.4 Обоснование выбора метода конечных элементов 42
2 Методика определения конечно-элементной модели гравитационного
поля Земли по измерениям 44
Постановка задачи и доказательство единственности её решения 44
Дискретизация области определения геопотенциала 54
Алгоритм поиска адреса конечного элемента 62
Преобразование координат и производных от возмущающего потенциала из общеземной в локальную и субпараметрическую системы координат 80
Аппроксимация возмущающего потенциала внутри субпараметрического конечного элемента алгебраическими полиномами 100
Определение базисных функций субпараметрического конечного конечного элемента 105
Вычисление коэффициентов и свободных членов линейных уравнений 115
Определение граничных условий в спутниковом шаровом слое и на земной поверхности 124
Построение траекторий движения ИСЗ 124
Вариационные уравнения для производных от векторов положения и скорости ИСЗ по узловым параметрам конечных элементов 139
Вариационные уравнения для изохронных производных 145
Аналитические формулы Суханова А. А. для вычисления изохронных производных 146
Параметрические уравнения поправок для дальностей и скорости изменения дальностей 147
Параметрические уравнения поправок для градиентометричес-ких измерений 150
Определение граничных условий на земной поверхности 160
3 Методика определения конечно-элементной модели гравитационного
поля Земли по заданному разложению геопотенциала в ряд по шаровым
функциям 181
Определение геопотенциала тяготения и его производных с помощью моделей разложения в ряд по шаровым функциям 181
Определение нормального потенциала силы тяжести и его производных 188
Вычисление возмущающего потенциала и его производных 198
Построение конечно-элементной модели 202
4 Описание компьютерных программ, алгоритмов и вычислительных
экспериментов 204
4.1 Программа построения спутниковых траекторий GENORX 204
Программа вычисления потенциала тяготения и его первых и вторых производных GRAVP 204
Программа вычисления высот квазигеоида GEOID 205
Система высот на основе модели EGM96 206
Программа вычисления уклонения отвесной линии UOL 208
Программа определения возмущающего потенциала и высот квазигеоида внутри нивелирного полигона методом конечных элементов МКЕ 210
Вычислительные эксперименты по определению возмущающего потенциала внутри нивелирного полигона методом конечных элементов 212
Вычислительный эксперимент по построению отдельного конечного элемента с использованием модели EGM96 214
Заключение 218
Список использованных источников 221
Приложение А Таблица - Координаты узлов икосаэдра 233
Приложение Б Таблица - Координаты узлов икосаэдра и додекаэдра 234
Приложение В Таблица - Координаты узлов икосаэдра и икосододекаэдра 237 Приложение Г Построение траектории ИСЗ GOCE по программе
GENORX 240
Приложение Д Высоты квазигеоида над эллипсоидом WGS-84 в
Новосибирской области 246
Приложение Е Определение возмущающего потенциала и высот квазигеоида внутри нивелирного полигона методом конечных
элементов по программе МКЕ 248
Приложение Ж Таблица - Значения возмущающего потенциала в узлах
конечного элемента 252
Введение к работе
Актуальность проблемы. Проблема построения глобальной высокоточной модели гравитационного поля Земли (ГПЗ) остается до сих пор нерешенной. Это связано с недостатком гравиметрической информации на большей части континентов, а также в приполярных областях. Достаточно плотной и качественной гравиметрической съёмкой покрыты лишь территории бывшего Советского Союза, США, Европы и Австралии. Кроме этого, серьёзной проблемой является согласование наземных гравиметрических сетей разных стран. Например, при сравнении гравитационных данных Советского Союза и Китая обнаружилось, что, в местах перекрытия, рассогласование этих данных огромно и достигает 70 мГал [1].
Альтиметрические измерения также не решают проблемы, так как ими не охвачены суша и приполярные зоны и по ним определяются не высоты геоида, а высоты морской поверхности (отличие достигает 1,5 м).
По траекторным спутниковым измерениям с наземных пунктов уверенно находятся лишь низкочастотные составляющие разложения ГПЗ в ряд по шаровым функциям. Следовательно, построить высокоточную модель гравитационного поля Земли только по наземным гравиметрическим, альтиметрическим и траекторным измерениям невозможно.
В настоящее время с появлением новых спутниковых градиентометриче-ских и спутник-спутниковых средств измерений, которые чувствительны к высокочастотной составляющей разложения ГПЗ в ряд по шаровым функциям, открылась возможность осуществить уточнение модели ГПЗ. Измерения уже выполняются, и будут в дальнейшем выполняться, глобально и плотно над всей поверхностью Земли на высотах от 240 до 500 км. Предлагаются и частично реализованы за рубежом три научных проекта, основанных на этих новых видах измерений [1-6]:
CHAMP (Challeging Minisatellite Payload) - совместный проект немецких, французских и американских ученых. Спутник запущен с космодрома Плесецк 15 июля 2000 г. Орбита круговая, высотой 450 км и с наклонением 87,27. Для определения ГПЗ используются измерения трёхосного бортового акселерометра и спутник-спутниковые измерения «высокий - низкий» (GPS-измерения). По шестнадцатисуточным данным уже построена модель OSU02A до пятидесятой степени и порядка. По точности эта модель, для волн длиной 800 км, на 10 % выше ранее полученных моделей, особенно в приполярных областях [7]. По трехмесячным акселерометрическим и спутник-спутниковым измерениям программы CHAMP построена также германо-французская модель длинноволнового гравитационного поля EIGEN-1S, которая по точности в два раза выше предыдущих моделей. Геоид, построенный по этой модели, характеризуется погрешностью не более 20 см. [8].
GRACE (Gravity Recovered and Climate Experiment) - проект, разработанный американскими учеными в Национальном управлении по аэронавтике и исследованию космического пространства (NASA). Пара искусственных спутников Земли (ИСЗ) проекта была запущена 17 марта 2002 г. с космодрома Плесецк. В апреле этого же года началось тестирование и калибровка научной аппаратуры ИСЗ. Полученные данные доказали работоспособность научной концепции проекта и подтвердили отличные характеристики измерительных устройств [9]. Спутники запущены на одну и туже близполярную орбиту высотой 500 км и на удалении друг от друга от 170 до 270 км. Расстояние между двумя ИСЗ измеряются с субмикронной точностью. На ИСЗ производятся GPS-измерения и акселерометрические измерения [10]. Погрешность измерения относительных ускорений 1 мкм/с, а погрешность GPS-измерений 1 см [11]. Предполагается, что, после обработки измерительной информации, модель статического гравитационного поля будет на несколько порядков точней существующих моделей [12].
GOCE (Gravity field and steady - state Ocean Circulation Explorer) - проект Европейского космического агентства [13-18]. Запуск планируется осуще-
ствить в ближайшем будущем с космодрома Плесецк. Вес спутника 1200 кг. Орбита ИСЗ близкруговая, солнечно-синхронная, с наклонением 96,5 и высотой от 240 до 270 км. Для точного определения орбиты на ИСЗ устанавливается высокоточный GPS-приемник. Основным инструментом для определения ГПЗ будет гравитационный градиентометр [19, 20], разработанный французской фирмой ONERF. Градиентометр представляет собой систему из шести трехосных электростатических акселерометров. Основой каждого акселерометра является пробная масса в 320 г, поддерживаемая в подвешенном состоянии внутри блока, вблизи его центра, с помощью электростатического поля, создаваемого электродами, размещенными на стенках блока. Перемещения этой массы ежесекундно измеряются. Результаты измерений фильтруются и преобразуются в гравитационные градиенты по восемнадцати различным направлениям. Точность измерения градиентов на уровне Ы0"3 Этвеш/^[Гц
[21]. Предполагается создать новую модель ГПЗ, по которой можно будет вычислять аномалии силы тяжести с погрешностью 1 мГал, а высоты геоида с погрешностью 1 см и с пространственным разрешением 70 км [18,22].
Для будущих космических проектов по изучению гравитационного поля Земли зарубежные ученые разрабатывают лазерный дальномер со стабилизированными длинами волн [23]. Он предназначен для измерения расстояний с
погрешностью 2,7-10" м/^Гц между двумя близкими ИСЗ, движущимися по
одной и той же орбите. По разности расстояний за определённые промежутки времени определяются вторые производные от геопотенциала по линии спутник-спутник.
В нашей стране также ведутся научно-исследовательские работы по использованию градиентометрических и спутник-спутниковых измерений для определения гравитационного поля Земли [13, 24]. Научно-производственное объединение прикладной механики предлагает включить в Федеральную космическую программу до 2010 года проект «Геопот» [25]. Он предназначен для
определения параметров ГПЗ и их изменений во времени. Масса ИСЗ 800 кг, высота орбиты 400 км, наклонение 90, срок службы - не менее семи лет.
Однако на высотах от 240 до 500 км, где выполняются спутник-спутниковые и спутниковые градиентометрические измерения, гравитационное поле имеет более гладкую структуру, чем у поверхности Земли. Поэтому, определяемые по спутниковым измерениям коэффициенты разложения ГПЗ в ряд по шаровым функциям будут отличаться от коэффициентов, которые были бы получены, если такого же типа измерения выполнить глобально на всей поверхности Земли. Отсюда возникает пока не решенная проблема продолжения вниз.
Частично эта проблема может быть решена, если поставить и решить краевую задачу с двумя граничными поверхностями. Одна граничная поверхность должна принадлежать спутниковому шаровому слою со спутниковыми градиентометрическими и/или спутник-спутниковыми измерениями. Другая граничная поверхность - земная поверхность, с максимально возможным набором наземных типов измерений, выполненных с достаточной точностью. Наземных измерений будет явно не хватать, поэтому в некоторых районах земной поверхности придется выполнять аппроксимацию и интерполирование функционалов от потенциала силы тяжести. Неизбежным является принятие каких-либо гипотез и допущений. Например, на море, где не выполнены гравиметрические съемки, считать морскую поверхность, вычисленную при помощи аль-тиметрических измерений, поверхностью геоида, с заданным значением потенциала силы тяжести Wo. На суше необходимо использовать нивелирные, гравиметрические, градиентометрические и GPS-измерения для вычисления потенциала силы тяжести и первых производных от него.
При построении моделей гравитационного поля Земли обычно решаются различные внешние краевые задачи теории потенциала в области, ограниченной земной поверхностью и бесконечной внешней границей. До недавнего времени целью исследований было совместное определение внутренней границы области, т. е. физической поверхности Земли и внешнего ГПЗ. В настоящее
время GPS-измерения и спутниковая альтиметрия позволяют, принципиально, новыми геометрическими методами решать первую часть задачи - определять геометрию земной поверхности. Актуальной проблемой остается вторая часть задачи - уточнение внешнего гравитационного поля Земли из совместной обработки наземных и новых спутник-спутниковых и спутниковых градиенто-метрических измерениий.
Гравитационное поле Земли можно описать глобально при помощи разложения в ряд по шаровым функциям, набором точечных масс или мультипо-лей, аномалиями силы тяжести и т. д. В рассмотренных выше космических проектах уже определяются и планируется дальнейшее уточнение коэффициентов Спт ,Snm разложения гравитационного поля в ряд по шаровым функциям. Недостаток такого разложения и других перечисленных способов представления гравитационного поля состоит в том, что для вычисления какой-либо характеристики поля в заданной точке на Земле или вне Земли с большой точностью требуется выполнить значительный объём вычислительной работы (просуммировать все члены ряда) и хранить большой объем информации. Кроме этого, при использовании математического аппарата гармонического анализа приходится относить выполненные на земной поверхности или в космосе измерения к сфере или решать проблему редуцирования измерений на сферу.
Избежать редуцирования и существенно сократить затраты на вычисления возможно, если разбить область определения ГПЗ на малые подобласти -конечные элементы и представлять в них гравитационный потенциал и функции от него локальными алгебраическими или тригонометрическими полиномами низкой степени. В этом случае придется сузить область определения геопотенциала и определять его в слое между земной поверхностью и спутниковым шаровым слоем. Сужение области - с одной стороны недостаток, но с другой стороны оно может повлиять на точность определения геопотенциала, потому что геопотенциальная функция, оказывается «зажатой» между двумя близкими поверхностями с известными значениями потенциала и его производных.
Преимущество таких локальных представлений гравитационного поля Земли по сравнению с глобальными, в задаче определения орбит спутников, показано в работах [27 - 29].
Если решать краевую задачу методом конечных элементов, то проще аппроксимировать сложную форму граничной поверхности Земли. Построенная путем решения краевой задачи конечно-элементная модель ГПЗ значительно удобней для практического применения, по сравнению с другими моделями. Из всего набора конечных элементов выбираются для использования только те ячейки, которые непосредственно попадают в район практической деятельности. Разработке такой конечно-элементной модели ГПЗ и посвящено данное исследование.
Какими оптимальными по форме и размерам должны быть конечные элементы? Какого вида, и какой степени должны быть полиномы? Как по координатам заданной точки определить в каком конечном элементе она находится? Каким образом подключить к решению краевой задачи разнородные наземные измерения и новые спутник-спутниковые и спутниковые градиентомет-рические измерения и определить значения потенциала в узлах конечных элементов? Как выполнить интерполирование и найти потенциал и его первые, и вторые производные внутри конечного элемента? Эти и другие актуальные вопросы являются предметом исследования диссертации.
Степень разработанности проблемы. Решению проблемы определения ГПЗ посвятили свои труды многие отечественные и зарубежные ученые: Молоденский М. С, Юркина М. И., Еремеев В. Ф., Бровар В. В., Шимбирев Б. П., Макаров Н. П., Пеллинен Л. П., Жонголович И. Д., Грушинский Н. П., Машимов М. М., Мещеряков Г. А., Марченко А. Н., Петровская М. С, Нейман Ю. М., Бывшев В. А., Остач О. М., Непоклонов В. Б., Плешаков Д. И., Бузук В. В., Вовк И. Г., Канушин В. Ф., Суздалев А. С, Панаев Г. А., Сурнин Ю. В., Кужелев С. В., Рапп Р., Чернинг С, Краруп Т., Шварц К., Руммель Р., Хейсканен В., Мориц Г., Кау-ла У. М., Бьерхаммар А., Бурша М. и др. Однако, решаемая проблема уточне-
ния моделей ГПЗ по новым спутник-спутниковым и спутниковым градиенто-метрическим измерениям является ещё недостаточно разработанной. Реализация проектов CHAMP и GRACE ещё не завершена, получены только первые предварительные результаты. Запуск ИСЗ программы GOCE ещё только планируется произвести и выполняются лишь модельные исследования. Трехмерные конечно-элементные модели ГПЗ, в области между земной поверхностью и спутниковым шаровым слоем, ранее не определялись.
Целью исследования являлась разработка методики определения конечно-элементной модели ГПЗ по разного вида спутниковым и наземным измерениям и по уже существующим разложениям геопотенциала в ряд по шаровым функциям.
Задачи исследования. Для достижения поставленной цели требовалось решить следующие задачи:
постановка краевой задачи;
доказательство единственности ее решения;
определение граничных значений производных геопотенциала в спутниковом шаровом слое;
определение граничных значений геопотенциала и его производных на земной поверхности;
дискретизация области определения геопотенциала;
замена интегральных уравнений системой линейных уравнений;
решение системы линейных уравнений для определения параметров модели ГПЗ;
вычисление геопотенциала и его производных по известным моделям разложения в ряд по шаровым функциям;
определение нормального потенциала силы тяжести и его производных;
- вычисление возмущающего потенциала в узлах конечных элементов.
Объект исследования - внешнее гравитационное поле Земли в области ограниченной земной поверхностью и спутниковым шаровым слоем.
Предмет исследования - методика построения конечно-элементной модели гравитационного поля Земли.
Теоретическая, методологическая и информационная база исследований. Рассмотрены различные численные методы решения задач математической физики. Из сопоставления методов, в качестве основы для решения проблемы уточнения ГПЗ, выбран метод конечных элементов. Экспериментальные исследования выполнены методом математического моделирования. Реальное поле заменено моделью гравитационного поля EGM96 до 360-й степени и порядка. Эта модель позволяет вычислить и считать «эталонными» значения потенциала и его производных в любой точке на и вне поверхности Земли, а также моделировать измерения, которые требуется обработать по предлагаемой методике. Достоинство метода математического моделирования состоит в том, что заранее известны параметры гравитационного поля, которые должны быть получены с помощью построенной конечно-элементной модели. Этот метод применяется для проверки правильности работы новых алгоритмов и программ, предназначенных для обработки реальных измерений. В диссертационной работе основные аспекты предлагаемой методики апробированы на модельных измерениях. Обработка реальных измерений возможна в будущем после того, как появится доступ к спутниковой и наземной измерительной информации.
Научная новизна работы:
внешний гравитационный геопотенциал определяется в области, ограниченной земной поверхностью и спутниковым шаровым слоем;
область определения геопотенциала дискретизируется шаровыми слоями и усечёнными треугольными пирамидами (конечными элементами) внутри слоев;
разработан алгоритм поиска номера конечного элемента, в котором расположена точка с заданными координатами;
получены формулы прямых и обратных преобразований координат точек и производных от возмущающего геопотенциала из общеземной системы
координат в локальную систему координат субпараметрического конечного элемента;
для определения граничных условий в спутниковом шаровом слое получены формулы связи спутник-спутниковых и спутниковых градиентометри-ческих измерений с возмущающим потенциалом в узлах конечных элементов;
для определения граничных условий на земной поверхности получены уравнения Лапласа в криволинейной геодезической системе координат и в системе плоских прямоугольных координат Гаусса - Крюгера и геодезической высоты;
разработана методика определения возмущающего потенциала внутри нивелирного полигона по спутниковым, гравиметрическим и градиентометри-ческим измерениям.
получены формулы для вычисления первых и вторых производных от нормального потенциала силы тяжести по прямоугольным пространственным координатам.
Теоретическая и практическая значимость работы. Предложенная в работе методика дает возможность определить удобную для практического использования трехмерную конечно-элементную модель ГПЗ. Разработанные для решения краевой задачи алгоритмы и компьютерные программы могут быть также использованы для вычисления орбит ИСЗ, высот квазигеоида и уклонений отвесных линий, для определения возмущающего потенциала внутри нивелирного полигона по гравиметрическим и градиентометрическим измерениям.
В первом разделе диссертации приведены общие сведения о потенциалах тяготения и силы тяжести, перечислены типы краевых задач, даны общие сведения о приближенных методах решения краевых задач. Рассмотрены: метод взвешенных невязок, метод коллокаций, метод конечных разностей, метод кол-локаций с подобластями, метод Бубнова - Галеркина, метод наименьших квадратов, метод конечных элементов, метод граничных интегральных уравнений, решение уравнения Лапласа методом разделения переменных, метод точечных масс. Показаны достоинства и недостатки перечисленных методов.
В результате анализа различных методов в качестве основы для исследований был выбран метод конечных элементов.
Во втором разделе рассмотрена методика определения конечно-элементной модели ГПЗ по измерениям, выполненным на земной поверхности и в спутниковом шаровом слое. Дана постановка задачи и доказана единственность её решения. Предложено дискретизировать область определения возмущающего потенциала шаровыми слоями и усеченными треугольными пирамидами (конечными элементами), построенными на основе многогранников. Получены формулы для вычисления базисных функций конечных элементов, а также коэффициентов и свободных членов линейных уравнений, из решения которых находятся значения возмущающего потенциала в узлах конечных элементов. Рассмотрены вопросы об определении граничных условий на внутренней поверхности спутникового шарового слоя и на земной поверхности.
В третьем разделе представлена методика определения конечно-элементной модели ГПЗ по заданному разложению геопотенциала в ряд по шаровым функциям. В основе методики лежат модифицированные автором формулы Пайнза С. для вычисления потенциала тяготения и его первых и вторых производных, а также формулы Молоденского М. С. для вычисления нормального потенциала. Они позволяют определять возмущающий потенциал в узлах конечных элементов и создавать конечно-элементные модели ГПЗ. Получены формулы для вычисления первых и вторых производных от нормального потенциала по направлению осей общеземной прямоугольной пространственной системы координат.
В четвертом разделе диссертации дано описание компьютерных программ, алгоритмов и вычислительных экспериментов.
Основные положения, выносимые на защиту следующие:
математическое описание конечно-элементной модели гравитационного поля Земли;
методика определения параметров конечно-элементной модели гравитационного поля Земли по измерениям из решения краевой задачи;
3) методика определения параметров конечно-элементной модели гравитационного поля Земли по заданному разложению геопотенциала в ряд по шаровым функциям.
