Специальный космофизический практикум:

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

В. М. Журавлев

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

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

Задача может быть рекомендована студентам физических факультетов университетов и других ВУЗов (в рамках общего практикума по механике).

1.   Гравитационное поле Земли и планет

1.1.   Гравитационное поле точечных масс

В основе теории тяготения Ньютона лежит закон всемирного притяжения, который для точечных масс величиной m₁ и m₂, находящихся на расстоянии r, может быть записан в виде

Здесь F₁₂ - сила, действующая на первое тело со стороны второго, r - радиус-вектор из первой точки ко второй, G = 6.67·10^-11 м³кг^-1с^-3 - ньютоновская гравитационная постоянная. Работа по удалению одного из тел в бесконечность при первоначальном расстоянии между точками r равна:

Эта функция определяет энергию взаимодействия точек и называется потенциальной энергией взаимодействия.

Величины g₂ = F₁₂/m₁ и φ₂ = - U₁₂/m₁ называются напряженностью гравитационного поля и его потенциалом, создаваемым второй материальной точкой в месте расположения первой точки. Эти величины не зависят от массы первой точки и характеризуют само поле вне зависимости от того, есть в данной области пространства материальная точка или нет. Заметим также, что напряженность и потенциал поля связаны между собой простым соотношением:

или в координатной записи

Кроме этих свойств, гравитационное поле обладает еще двумя свойствами. Одно из них - принцип линейной суперпозиции, а второе - отсутствие экранирования. Первое можно сформулировать следующим образом. Пусть имеются n материальных точек с массами m₁, m₂,:,m_n, расположение которых в пространстве определяется набором радиус-векторов r₁, r₂,:,r_n. Тогда гравитационное поле, создаваемое этими материальными точками в точке с радиусом-вектором r₀ характеризуется напряженностью и потенциалом, которые могут быть получены в виде линейной суммы отдельных напряженностей и потенциалов

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

Рассмотрим предельный переход к сплошной среде. В этом случае каждая точка среды с координатами r' = (x', y', z') может рассматриваться как материальная точка, имеющая массу dm(r') = ρ(r')dV', где ρ(r') - плотность распределения массы в данной точке, а dV' - элементарный объем среды в этой точке. Тогда в точке пространства с координатами r = (x, y, z) данный элемент среды создает гравитационное поле со следующими параметрами:

Используя принцип суперпозиции, получаем, что в точке r суммарное поле характеризуется следующими величинами:

		(1)

		(2)

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

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

1.2.   Гравитационное поле тел общего вида

Общий вид потенциала, создаваемого телом произвольной формы дается формулой (2). В эту формулу входит множитель

где ψ - угол между векторами r и r'. Правую часть последнего соотношения можно разложить в ряд по полиномам Лежандра:

		(3)

где

Кроме этого введем присоединенные полиномы Лежандра, которые определяются следующими соотношениями:

В формулу (3) входит угол ψ, который выражается через сферические углы θ и φ системы координат. Для преобразования (3) к виду, в который входят θ и φ, воспользуемся теоремой сложения для полиномов Лежандра. Согласно этой теореме

Здесь θ и φ - сферические углы вектора r, направленного в точку наблюдения, а θ' и φ' - сферические углы вектора r', направленного в точку расположения текущего источника. Подставляя последнее выражение в формулу (2) приходим к соотношению:

		(4)

1.3.   Гравитационное поле планет

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

Если тело слабо деформировано, то ряд (4) будет содержать малые поправки к ньютоновскому полю. Планеты и вращающиеся одиночные звезды обладают осевой симметрией: ρ = ρ(r, θ). В этом случае секториальные и тессеральные составляющие гравитационного поля равны нулю:

В результате ряд (4) упрощается и приобретает вид

		(5)

В случае зеркальной симметрии тела относительно экваториальной плоскости нечетные составляющие в ряде (5) исчезают. Последнее условие выполняется для одиночных вращающихся звезд. В этом случае ряд (5) можно записать так:

		(6)

Таблица 1. Данные о зональных гармониках планет солнечной системы [2]

		Jn·106
Небесное тело	GM	n
	м3/с2	1	2	3	4	5
Солнце	1.327·1020
Меркурий	2.168·1013
Венера	3.249·1014	-20
Земля	3.987·1014	-1082	2.54	1.58	0.22	-0.59
Луна	4.903·1012	-206	-37.7	33.3	-5.5
Марс	4.298·1013	-1995	29	-9.5	5	0.2
Юпитер	1.267·1017	-14750		580		5
Сатурн	3.793·1016	-16670		1040		11
Уран	5.803·1015	-12000
Нептун	7.026·1015	-3900

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

		(7)

2.   Движение точечной частицы в поле тяготения Земли

2.1.   Кеплеровские орбиты

Уравнения Ньютона движения точечной частицы массой m в поле тяготения имеют вид

		(8)

где φ - гравитационный потенциал поля тяготения. Для движения в поле сферически симметричного тела потенциал имеет вид:

В сферическом поле сохраняется момент импульса:

В силу этого удобно выбрать систему координат таким образом, что бы плоскость движения спутника (плоскость орбиты) была бы ортогональна постоянному вектору L. Направим ось z вдоль вектора L. В этом случае L = (0, 0, L). Переходя теперь к полярным координатам в плоскости орбиты получаем

		(9)

Рис. 1. Элементы эллиптической орбиты.

Вторым интегралом движения является интеграл энергии, который в выбранной полярной системе координат имеет вид

		(10)

Здесь v_r - радиальная скорость спутника, E₀ - полная энергия. Учитывая (9), приходим к следующему уравнению:

		(11)

Делая замену переменных r(t) = r(φ(t)) получаем:

Далее делаем замену ξ(φ) = r ^-1(φ) и в результате получаем:

или

Этот интеграл совпадает по форме с интегралом гармонического осциллятора единичной массы и единичной частоты и точкой равновесия равной

Поэтому можем сразу записать решение:

Обычно это решение записывается в следующей форме

		(12)

где p - фокальный параметр, e - эксцентриситет, угол ν = φ - φ₀, отсчитываемый от направления на перигей, называется истинной аномалией. Перигеем называется точка минимального расстояния орбиты от центра поля, апогеем - точка максимального удаления, если она есть. Угол поворота относительно центра эллипса, отсчитываемый от перигея, называется эксцентрической аномалией E. В апогее и перигее радиальная компонента скорости равна нулю. Поэтому из интеграла энергии находим:

Отсюда

С другой стороны из (12) имеем

		(13)

Отсюда

		(14)

Движение спутника по орбите определяется с помощью уравнения Кеплера, которое имеет следующий вид:

		(15)

Здесь E - эксцентрическая аномалия, величина n(t - t₀) + ₀ называется средней аномалией, а величина ₀ - средней аномалией в эпоху (t₀), - среднее движение. Геометрический смысл эксцентрической аномалии пояснен на рис. 1. Эксцентрическая аномалия связана с естественной аномалией (угол ν на рис. 1) формулой

		(16)

Из этих двух уравнений определяется истинная аномалия, соответствующая определенному моменту времени t и, следовательно, положение спутника на орбите, например, относительно перигея.

Декартовы координаты спутника в геоцентрической системе координат будут иметь при этом следующий вид:

		(17)

		(18)

		(19)

Здесь u = ν + ω.

Подробности можно найти в [1]. Полный набор формул смотрите в [5] (стр. 171-180).

2.2.   Положение орбиты в пространстве

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

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

Рис. 2. Положение орбиты в пространстве.

Угол в плоскости орбиты, отсчитываемый от восходящего узла до перигея называется аргументом перигея и обозначается ω, угол между плоскостью орбиты и экваториальной плоскостью - наклонением орбиты и обозначается здесь через i. Так же в параметры орбиты входит фокальный параметр p и эксцентриситет e. Эти параметры приведены на рис. 2.

2.3.   Движение в сфероидальном поле. Модель прецессии

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

		(20)

где m - масса спутника, и использовано тождество cosθ = z/r. Уравнения Ньютона движения спутника в таком поле имеют вид

Умножая это уравнение слева векторно на радиус-вектор r = xe_x + ye_y + ze_z, где e_x, e_y, e_z - орты декартовой системы координат, получаем уравнение для момента импульса спутника L = [r × p] в следующем виде:

		(21)

Если сила, действующая на спутник центральна, то момент сил, стоящий в правой части (21) равен нулю. Действительно, для центральной силы выполняется условие:

где Q(r) - некоторая скалярная функция длины радиус вектора r = |r|. Для сферической планеты гравитационное поле центрально и момент сил равен нулю. Это приводит к закону сохранения момента импульса. Если форма Земли отлична от сферической, то гравитационное поле не является центральным, и момент импульса не сохраняется. Вычислим момент сил для потенциальной энергии (20). Дифференцируя потенциальную энергию (20) получаем следующее выражение для силы, действующей на спутник:

		(22)

где

а вектор e_z - единичный вектор в направлении оси z, совпадающей с направление оси вращения Земли. Отсюда следует, что в таком поле момент сил, действующих на спутник равен

Отсюда находим

Учитывая, что

получаем для компонент момента импульса систему уравнений

		(23)

Из этой системы видно, что проекция момента импульса спутника на ось вращения Земли сохраняется: L_z = const.

Уравнения (23) вместе с законом сохранения энергии теперь полностью определяют движение спутника в рассматриваемом поле тяготения. Поскольку отклонение от сферичности у Земли мало (см. таб. 1), можно считать, что за один оборот спутника вокруг Земли существенных изменений параметров орбиты не происходит. Такие отклонения должны проявляться лишь через значительный промежуток времени. Поэтому можно предположить, что закон орбитального движения спутника остается по форме таким же, но параметры орбиты медленно меняются со временем. Это означает, что для вычисления декартовых координат спутника можно использовать те же соотношения (17), но теперь следует считать Ω = Ω(t), ω = ω(t), i = i(t), ε = ε(t), a = a(t). Всю совокупность уравнений для этих параметров орбиты, которые называются оскулирующими элементами, можно получить теперь из уравнений (23) и закона сохранения энергии. Однако эти вычисления громоздки, и мы их полностью приводить не будем. Приведем лишь вывод уравнения для векового изменения долготы восходящего узла Ω в предположении, что орбита имеет малый эксцентриситет, т. е. практически не отличается от круговой.

Подставим выражения для координат x, y, z из (17) в первые два уравнения (23). В результате получим

Отсюда получаем

		(24)

		(25)

В этих уравнениях в правой части имеются слагаемые, которые меняются быстро за один оборот спутника вокруг Земли и те которые меняются медленно. Быстро меняются слагаемые, содержащие истинную аномалию ν, которая за один оборот меняется от 0 до 2π. Медленно меняются слагаемые, не содержащие истинной аномалии. Из этого анализа следует, что выражения для компонент момента импульса можно представить аналогичным образом в виде слагаемых, меняющихся быстро и медленно, т. е.

где функции L_x⁽⁰⁾, L_y⁽⁰⁾ меняются медленно. Уравнения для этих функций получаются прямой подстановкой последних выражений для компонент момента импульса в уравнения (24) и (25) и приравниванием медленно меняющихся частей уравнения. В результате имеем