Размещено на http://www.allbest.ru/

Размещено на http://www.allbest.ru/

Введение

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

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

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

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

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

Все вычисления производятся в среде MATLAB. Для анализа временного ряда колебаний J₂ используются Фурье анализ и сингулярный спектральный анализ (ССА).

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

В данной работе используется временной ряд колебаний коэффициента J₂ для Земли, измеренный с 1976 года до марта 2016 года. Данные о планетах Солнечной системы взяты с сайта NASA, последнее обновление этих данных состоялось 28 января 2016 года. Коэффициенты, не представленные на этом сайте, были взяты из последних печатных работ, в которых они вычислялись на основе наиболее современных небесномеханических моделей.

Вычисления по формулам Клеро и Радо-Дарвина проводил в своих книгах В.Н. Жарков [1]. Колебания коэффициента гравитационного поля J₂ исследовались в работе М. Ченга [2], где использовался метод вейвлет-анализа. Также J₂ исследовалось с помощью Фурье анализа [3]. Метод сингулярного спектрального анализа применялся для изучения других процессов, протекающих на Земле, к примеру, изменения содержания озона в атмосфере планеты [4] или колебаний глобальной температуры [5]. В данном исследовании проводится сравнение найденных компонентов колебаний коэффициента геопотенциала J₂ с другими процессами, происходящими на Земле. В работах [5,13] была предположена связь между колебаниями скорости вращения Земли и глобальными изменениями климата; в представленном нами исследовании с этими процессами сравниваются и колебания гравитационного поля.

Новизна данного исследования заключается в первую очередь том, что для анализа колебаний гравитационного поля применяется сингулярный спектральный анализ. Данный метод может способствовать нахождению новых, ранее необнаруженных трендов и периодичностей в исследуемом процессе. Также новым в этой работе является то, что в ней используются наиболее полные и наиболее современные данные о гравитационных полях планет Солнечной системы. Что касается временного ряда для колебаний коэффициента J₂, то в большинстве предыдущих исследований используются данные, полученные начиная с конца 80-х годов прошлого века, тогда как данное исследование использует измерения, производившиеся с 1976 года до весны 2016, что является наиболее полными измерениями на сегодня.

В результате данного исследования проверено соответствие фигур планет Солнечной системы гидростатически равновесным моделям. Разделены и отфильтрованы от шумов компоненты изменчивости коэффициента J₂ для Земли.

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

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

1.1 Основные сведения о гравитационном поле и фигуре планеты

Сжатие планеты

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

где a - экваториальный радиус, b - полярный радиус планеты. Правда, полученное им число б=1/230 (?0.00435) было ещё весьма неточным.

Заключение Ньютона о сжатии Земли оспаривалось многими учёными, в числе которых был знаменитый французский астроном Кассини Ж.Д. Для проверки того, сжата Земля у полюсов, или вытянута, в середине XVIII в. академией наук Франции были организованы экспедиции для выполнения градусных измерений на различных широтах нашей планеты. В результате проделанных измерений было подтверждено, что фигура Земли представляет собой сплюснутый сфероид, полярная ось которого примерно на 20 км меньше экваториальной. Точка зрения Ньютона о сплюснутости фигуры Земли получила экспериментальное подтверждение, благодаря чему она была признана другими учёными.

Моменты инерции планеты

Тензор инерции твёрдого тела задаётся матрицей:

- главные (диагональные) моменты инерции,

- центробежные моменты (произведения) инерции.

В системе координат главных осей инерции он приобретает диагональный вид:

По главной оси расположены осевой С и экваториальные А и B моменты инерции.

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

Пренебрегая трехосностью Земли, можно считать А=B.

Cредний момент инерции при этом даётся выражением

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

Важным дополнительным параметром для определения соотношения между C и А, помимо J₂, является постоянная прецессии:



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

Распределение плотности в недрах планеты существенно влияет на средний момент инерции I и, наоборот, значение I, определённое экспериментально, контролирует распределение плотности при модельных расчётах.

Введём ещё одну важную характеристику - безразмерный момент инерции

В случае планеты постоянной равномерно распределённой плотности её безразмерный момент инерции I^* равен 0.4. Легко убедиться путём непосредственных численных расчётов, что при росте плотности в недрах планеты от периферии к центру величина I? будет принимать значение, меньшее 0.4. Наоборот, если в планете происходит уменьшение плотности c глубиной, то значение I? будет превосходить предельное значение, равное 0.4 (полых планет пока не обнаружено).

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

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

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

Гравитационный потенциал и коэффициент J2

Потенциалом принято называть работу, которую нужно совершить, чтобы переместить материальную точку единичной массы из заданной точки в бесконечно удалённую. Пусть F есть вектор силы (записан в столбец), приложенной к материальной точке, r - радиус-вектор этой точки. Тогда потенциалом в данной точке будет величина:

Выражение под интегралом в данной формуле является полным дифференциалом силовой функции, т.е.

поэтому:

т.к. силовую функцию в бесконечно удалённой точке можно приравнять нулю (V(?) = 0).

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

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

Выберем систему координат так, чтобы одна из материальных точек оказалась в начале координат этой системы. Тогда другая материальная точка будет иметь радиус-вектор r.

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

Где - гравитационная постоянная. Абсолютная величина вектора F равна:

Соответственно, гравитационный потенциал точки (силовая функция) равен:

где - расстояние между притягивающимися точками (скалярная величина).

Исходя из принципа суперпозиции, гравитационный потенциал точек равен сумме их гравитационных потенциалов:

Пусть точек бесконечно много, а массы их бесконечно малы. Тогда:

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

Пусть - координаты точки P, а о, з, ж - координаты текущей точки с массой dm; тогда гравитационный потенциал можно переписать в виде следующей формулы:

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

где r - расстояние от центра сферы, GM - планетоцентрическая постоянная.

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

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

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

Нечетные члены отсутствуют, поскольку поле такой планеты симметрично относительно экватора (теорема Лихтенштейна).

Ограничиваясь первым членом суммы, получаем

- гравитационный момент, a - экваториальный радиус Земли,

- второй полином Лежандра, A и С - главные моменты инерции, и - полярный угол, равный дополнению широты до , т.е. .

До запусков ИСЗ за счёт наземных измерений удалось определить первый поправочный член J₂ к основной части гравитационного потенциала Земли. В спутниковую эпоху точность определения существенно возросла.

1.2 Характеристики планет Солнечной системы

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

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

За последние годы благодаря космическим аппаратам, совершавшим пролёты мимо планет Солнечной системы, были получены уточнённые данные об особенностях их гравитационного поля и фигуры. Так, к примеру, автоматическая межпланетная станция «MESSENGER», вышедшая в 2011 году на орбиту Меркурия, позволила уточнить данные о гравитационном поле этой планеты, которые имелись у учёных на тот момент.

Таблица 1. Данные о коэффициенте гравитационного потенциала J2, безразмерном моменте инерции I* и сжатии планет Солнечной системы и Луны [6], [7], [8], [9]

Планета	,	Безразмерный момент инерции	Сжатие б
Меркурий	50.3	0.35	0
Венера	4.458	0.33	0
Земля	1082.63	0.3308	0.00335
Луна	202.7	0.394	0.0012
Марс	1960.45	0.366	0.00589
Юпитер	14736	0.254	0.06487
Сатурн	16298	0.210	0.09796
Уран	3343.43	0.225	0.02293
Нептун	3411	0.290	0.01708
Плутон	90.1	0.310	0

Значения сжатия планеты, равное 0, означает, что на современном уровне точности его определения планету можно считать сферической.

1.3 Расчёт сжатия б по формуле Клеро и сравнение с наблюдаемыми значениями

Теория фигуры планет во многом ведёт своё начало от работы французского математика А.К. Клеро «Теория фигуры Земли, основанная на началах гидростатики», опубликованной в 1743 г. Основываясь на законе всемирного тяготения, Клеро показал, что ускорение силы тяжести на поверхности Земли как функция широты изменяется по закону:

где ? - широта места, - ускорение силы тяжести на экваторе, (- отношение центробежной силы к силе тяжести на экваторе, - угловая скорость вращения Земли, a - её экваториальный радиус).

Теорема Клеро позволяет определить сжатие б независимо от определения геометрических элементов путём градусных измерений. Согласно теории Клеро, чтобы определить сжатие планеты б, достаточно определить гравитационной потенциал на её поверхности.

Соответственно, сжатие земного сфероида б простым образом связано с J₂, угловой скоростью вращения Земли щ, полной массой M и экваториальным радиусом a:

где - период обращения планеты.

Для проверки выполнения теории Клеро нами были проведены расчёты по формуле (1) с использованием известных на сегодня значений параметров планет. Для расчётов потребовались величины из таблицы 2:

Таблица 2. Экваториальный радиус, планетоцентрическая постоянная и наблюдаемый период обращения планет Солнечной системы и Луны

Планета	Экваториальный радиус , км	Планетоцентрическая постоянная ,	Период обращения планеты , часы
Меркурий	2439.7	0.022032	1407.6
Венера	6051.8	0.32486	5832.5
Земля	6378.137	0.39860	23.9
Луна	1738.1	0.00490	655.7
Марс	3396.2	0.042828	24.6
Юпитер	71492	126.687	9.9
Сатурн	60268	37.931	10.7
Уран	25559	5.7940	17.2
Нептун	24764	6.8351	16.1
Плутон	1187	0.000870	153.3

Полученные значения таковы:

Таблица 3. Вычисленное значение q, а также сравнение расчётной и наблюдаемой величины сжатия планет Солнечной системы и Луны

Планета		Расчётное б	Наблюдаемоеб
Меркурий	1.013 * 10-6	0.00076	0
Венера	6.109 * 10-8	0.000007	0
Земля	0.00347	0.00336	0.00335
Луна	7.592286 * 10-6	0.000308	0.0012
Марс	0.004604	0.005243	0.00589
Юпитер	0.0896447	0.066926	0.06487
Сатурн	0.153551	0.101223	0.09796
Уран	0.02967228	0.019851	0.02293
Нептун	0.0261108	0.018172	0.01708
Плутон	0.000249174	0.00026	0

Видно, что для всех планет Солнечной системы результат вычислений практически совпал с наблюдаемым значением сжатия. Единственное отличие есть только в случае с Луной.

1.4 Расчёт периода обращения по формуле Дарвина-Радо и сравнение с наблюдаемыми значениями

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

Исследования показали, что потенциал притяжения гидростатически равновесной планеты содержит лишь чётные зональные гармоники:

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

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

Формула Радо-Дарвина относительно периода:

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

Планета	Средняя плотность,	Полученное, часы	Нынешнее , часы
Меркурий	5427	70.5	1407.6
Венера	5243	220.8	5832.5
Земля	5514	13.9	23.9
Луна	3340	52.7	655.7
Марс	3933	14.1	24.6
Юпитер	1326	4.7	9.9
Сатурн	687	2.6	10.7
Уран	1271	6.7	17.2
Нептун	1638	11.5	16.1
Плутон	2095	70.4	153.3

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

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

1.5 Современные наблюдения за J2. SLR-метод

Лазерная локация спутников (SLR, Satellite Laser Ranging) - это широко применяемая в геодезии методика определения орбит. На основе изменений расстояний до спутника с установленными на его борту лазерными отражателями со множества станций на Земле определяются координаты спутника. Наносекундные импульсы лазера позволяют осуществить за один сеанс миллиарды измерений дальности, и на основе обработки таких измерений точность определения координат спутника достигает нескольких сантиметров и даже миллиметров.

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

Таким образом, чем точнее были определены координаты спутника, тем точнее можно определить параметры гравитационного поля, в том числе коэффициент гравитационного потенциала J₂.

На данный момент на орбите Земли работает несколько спутников, осуществляющих лазерную локацию планеты. На базе их наблюдений, находящихся в открытом доступе [10], [11], был составлен временной ряд для отклонений коэффициента J₂ от базового значения (рис. 1):

Рисунок 1. Вариации коэффициента гравитационного потенциала J2, полученные методом лазерной локации за период с 1976 по 2016 год

Миссия GRACE

Миссия GRACE (Gravity Recovery And Climate Experiment) специально запущенная в 2003 г. для измерений гравитационного поля Земли [24]. Однако пара спутников-близнецов GRACE определяет коэффициент C20 менее точно, чем спутники системы SLR, поэтому при обработке результатов миссии GRACE полученные значения C₂₀ (пропорционален - J2) заменяются значениями, полученными с помощью метода лазерной локации. На рис 2. представлены для сравнения коэффициенты Стокса C₂₀ разложения гравитационного поля, полученные по данным GRACE и SLR.



Рисунок 2. Сравнение коэффициента С20, полученного по данным GRACE и лазерной локации спутников SLR

2. Анализ временного ряда коэффициента J2 для Земли

2.1 Фурье-анализ

Описание метода

Для непрерывного преобразования Фурье используется выражение:

где - циклическая частота.

Для каждой заданной частоты результат Фурье-преобразования определяет меру присутствия гармоники частоты в сигнале:

Квадрат модуля преобразования Фурье называют спектральной плотностью, модуль - амплитудным спектром.

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

Автоковариационная функция (АКФ) (смещённая оценка) определяется по формуле

где fc - отцентрированные (вычтено среднее) значения временного ряда, заданные в моменты времени t_i, всего N значений, - задержка.

Для случайных процессов доказана теорема Винера-Хинчина об их представимости спектральной плотностью мощности СПМ. Для получения СПМ необходимо выполнить преобразование Фурье АКФ.

Исследование временного ряда колебаний J2

С помощью функций среды MATLAB был проведён анализ SLR ряда для J2.

На первом этапе была вычислена автоковариационная функция

Рис. 3. Автоковариационная функция исходного ряда для положительных задержек (единица - в месяцах)

Рисунок 4. Спектральная плотность мощности, полученная Фурье-преобразованием АКФ

Затем был проведён непосредственный Фурье-анализ исходного ряда с помощью БПФ (функция fft).

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

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

Колебательные процессы заметны на более коротких периодах:

Рисунок 6. Амплитудный спектр исходного ряда с выделенными пиками на некоторых из периодов

Самые большие из заметных периодов - примерно 10 и 20 лет.

Рассмотрим спектр под ещё большим увеличением:



Рисунок 7. Амплитудный спектр исходного ряда с выделенными пиками на некоторых из периодов

Здесь заметен ряд слабых компонент с периодами примерно 5.7, 4.5, 3.7, 3, 2.1 года. Однако не все из них достоверны из-за квази-случайной природы сигнала. Наиболее значительная по амплитуде компонента и совершенно достоверная компонента имеет период 1 год.

Рисунок 8. Короткие периоды амплитудного спектра исходного ряда с выделенными пиками на некоторых из периодов

На коротких периодах полученные амплитуды колебаний незначительны, выделяется только колебание с периодом 0,5 года.

2.2 Сингулярный спектральный анализ временного ряда

Описание метода

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

Алгоритм ССА

Пусть сигнал представлен временным рядом f(t_k), содержащим Nотсчётов. ССА выполняется в 4 этапа, первые два служат для разложения, другие два для восстановления.

Вложение

Выбирается число L (длина окна, параметр задержки, лаг) и формируется траекторная матрица, столбцы которой {xi} представляют собой последовательно выбранные элементы из временного рядя вектора размерности L:

В траекторное матрице Lстрок и K = N - L + 1 столбцов. В ней по боковым диагоналям, т.е. где стоят одинаковые элементы.

Сингулярное разложение траекторной матрицы:

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

Размерности матриц U = (LxL), V - (KxK). Из разложения квадратной матрицы

видно, что -её собственные числа, а -собственные вектора, являющиеся строками унитарной матрицы .

Группировка

Если матрица поворота Sсодержит dположительных сингулярных чисел , то ранг Х равен d и каждому сингулярному числу можно сопоставить тройку и компоненту разложения

Сгруппировав специальным образом тройки, полученные при разложении, можно представить исходную матрицу в виде

где - группа, содержащая компоненты с индексами . Таким образом, множество компонент {1, 2, 3, …, d} разбивается на mнепересекающихся подмножеств, в которых осуществляется группировка.

Генкелизация

Если сигнал удачно сгруппирован, каждая компонента будет представлять собой некоторую аддитивную компоненту f. Извлечём компоненту сигнала , содержащуюся в матрице , усреднением элементов вдоль побочных диагоналей i + j = k + 2. Обозначив компоненты как , и числа L* = min (L, K), K* = max (L, K), получим отсчёты по формулам

Применяя эту операцию для каждой , получаем представление исходного ряда в виде суммы m рядов , являющиеся главными компонентами (ГК) сигнала.

Исследование временного ряда колебаний J2

Ниже представлены результаты сингулярного спектрального анализа временного ряда ДJ₂ для параметра L = 18 лет (при месячном шаге данных L=216). Сравнение производилось также для L = 5, 10, 12 и 20 лет и был сделан вывод об оптимальности L=18 лет для наиболее чистого разделения компонент.

В анализе были использованы первые 10 сингулярных чисел. Далее продемонстрированы графики 10 полученных для них компонент. Для удобства восприятия мы разделили все компоненты на два отдельных графика (рис. 8, 9.):



Рисунок 9. Графики компонент, соответствующих сингулярным числам с 1 до 5

Рисунок 10. Графики компонент, соответствующие спектральным числам с 6 по 10

Для некоторых компонент (к примеру, СЧ1 и СЧ6) заметны краевые эффекты (резкие скачки по краям), что является особенностью используемого метода.

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

· Объединить СЧ1 и СЧ2, так как они полностью совпадают по периоду и практически совпадают по амплитуде;

· Объединить СЧ3 и СЧ4, так как они представляют собой практически идентичные линии тренда;

· Объединить СЧ6 и СЧ7, так как они полностью совпадают по периоду и практически совпадают по амплитуде;

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

После произведения указанных действий получен следующий результат:

Рисунок 11. Графики итоговых главных компонент, полученных после удаления шума и слияния некоторых спектральных чисел

Если просуммировать графики на рисунке 6 сравнить их сумму с графиком исходного ряда, получится следующий график:

Рисунок 12. График исходного временного ряда после удаления шума

Видно, что ряд, включающий первые 9 СЧ практически полностью выбирает дисперсию (изменчивость) исходного ряда и воспроизводит его с высокой степенью приближения. СЧ с номерами больше 9 го считаются в нашем исследовании шумовыми. График для объединённых 9 компонент является более сглаженным, на нём нет резких скачков. Это связано с тем, что из исходного ряда, отбрасыванием СЧ с номерами >9, был удалён шум. ССА является хорошим методом фильтрации шумов.

2.3 Анализ полученных главных компонент

Рассмотрим подробнее найденные компоненты ряда.

ГК1 (СЧ1+2)

График компоненты ГК1:

Рисунок 13. График компоненты, полученной после слияния спектральных чисел 1 и 2

С помощью Фурье-анализа был получен спектр этой компоненты (рис. 13). Представленный спектр отображён на периодов, благодаря чему можно увидеть период найденного колебания:

Рисунок 14. Квадрат и модуль амплитуды, полученные быстрым преобразованием Фурье ГК1

На графике явно видно, что период этой компоненты равен одному году. Амплитуда

На протяжении года минимума эта компонента достигает в январе-феврале, максимума в июле-августе.



Рисунок 15. Колебание ГК1 за 1 год

Изменение сжатия в годовом цикле связано с сезонными перераспределениями атмосферных и океанических масс. Снег и вода на континентальной поверхности также оказывает заметное воздействие на сжатие планеты. [12]

ГК2 (СЧ 3+4)

График компоненты ГК2:

Рисунок 16. График компоненты, полученной после слияния спектральных чисел 3 и 4

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

Если предположить, что в 70-е годы данная компонента находилась в максимуме, то её период может быть оценён в 60-70 лет.

Интерес представляет тот факт, что такого рода 60-летняя изменчивость присутствует также в огибающей Чандлеровского движения полюса Земли и длительности суток LOD [5,13].

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

ГК3 (СЧ 5)

График компоненты ГК3:

Рисунок 17. График компоненты, соответствующей спектральному числу 5

Период данной компоненты достаточно нестабилен. Использование Фурье-преобразования даёт следующий результат:



Рисунок 18. Квадрат и модуль амплитуды, полученные быстрым преобразованием Фурье ГК3

Период этой компоненты около 20-ти лет. Однако он выделяется недостаточно точно в связи с его малой амплитудой.

Найденная компонента гравитационного может соответствовать приливной компоненте с периодом в 18.6 года, найденной в [2]. Амплитуда месячных колебаний склонения Луны изменяется с периодом 18,61 г. от 29° до 18°, из-за прецессии лунной орбиты.

ГК4 (СЧ 6+7)

График компоненты ГК4:

Рисунок 19. График компоненты, полученной после слияния спектральных чисел 6 и 7

Периодический спектр этой компоненты, полученный с помощью Фурье-преобразования:

Рисунок 20. Квадрат и модуль амплитуды, полученные быстрым преобразованием Фурье ГК4

На графике явно видно, что эта компонента совершает 2 колебания в год, т.е. её период равен 0.5 года. Амплитуда порядка 5

На протяжении года минимума эта компонента достигает в мае и октябре, максимума в феврале и июле.

Рисунок 21. Колебание ГК4 за 1 год

Полугодовая компонента также связана с атмосферным и океаническим перераспределениями массы.

ГК5 (СЧ 8)

График компоненты ГК5:

Рисунок 22. График компоненты, соответствующей спектральному числу 8

Фурье-анализ этой компоненты позволяет определить ее период в 10 лет:

Рисунок 23. Квадрат и модуль амплитуды, полученные быстрым преобразованием Фурье ГК5

Амплитуда этого колебания приблизительно равна

Минимумы этой компоненты наблюдались в 1985, 1995 и 2005 годах, максимумы в 1990, 2000 и 2010. По краям графика из-за краевых эффектов отображение может быть недостаточно достоверным, однако периоды 1976-1985 и 2010-2016 также практически точно соответствуют обнаруженной закономерности.

В [15] обсуждается взаимосвязь изменения температуры за Земле и активность Солнца. Как известно, солнечная активность меняется в 11-ти летнем цикле, однако минимумы и максимумы компоненты J₂, представленной на рис. 21, не совсем синхронны с циклом солнечной активности.

Найденная компонента скорее связана с полупериодом прецессии орбиты Луны, равному 9.3 года.

ГК6 (СЧ 9)

График компоненты ГК6 отдельно от остальных:

Рисунок 24. График компоненты, соответствующей спектральному числу 9

Фурье-анализ этой компоненты позволяет выделить период в 8 лет:

Рисунок 25. Квадрат амплитуды, полученный быстрым преобразованием Фурье ГК6

Амплитуда этого колебания приблизительно равна

Минимумы этой компоненты наблюдались в 1987, 1995, 2004 и 2011 годах, максимумы в 1991, 1999 и 2008. По краям графика из-за краевых эффектов отображение может быть недостаточно достоверным, однако периоды 1976-1987 и 2011-2016 также практически точно соответствуют обнаруженной закономерности.

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

2.5 О влиянии изменений J2 на частоту Чандлеровского колебания

Известно, что изменения в J2 должны приводить к изменениям Чандлеровской частоты fc нутации земной оси (период~433 суток) согласно формуле:



где = 7.292 115 *10^-5 рад/c - средняя угловая скорость вращения Земли, С, А - моменты инерции R - экваториальный радиус G - гравитационная константа, k - число Лява.

Считая, что меняется лишь величина C и связав эти изменения с выражением для J₂, получаем для Чандлеровского периода:

Для T - 433 дня секунд

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

Заключение

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

· Был произведён расчёт сжатия планет Солнечной системы по формуле Клеро. Получено, что при вычислении по этой формуле получается результат, близкий к наблюдаемым параметрам. Однако, было обнаружено некоторое отклонение при расчётах для Луны - наблюдаемое сжатие отличается от расчётного почти в 4 раза, тогда как в остальных случаях соотношение очень близко к 1.

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

· Был произведён анализ SLR-временного ряда вариаций коэффициента гравитационного потенциала J₂ c 1976 по 2016 год методами Фурье анализа и сингулярного спектрального анализа.

· Были проанализированы компоненты, выделенные при помощи сингулярного спектрального анализа. Для получения периода этих компонент использовался Фурье анализ:

Таблица 5. Итоги сингулярного спектрального анализа временного ряда коэффициента J2

ГК №	Примерный период, г	Амплитуда	Минимумы	Максимумы
1	1		январь-февраль ежегодно	Июль-август ежегодно
2	60		2005	1970-е
3	20	-	1985, 1995 и 2006	1980, 2000, 2010
4	0.5		Май и октябрь ежегодно	февраль и июль ежегодно
5	10		1985, 1995 и 2005	1990, 2000 и 2010
6	8		1987, 1995, 2004 и 2011	1991, 1999 и 2008

· Найденные периодичности были сопоставлены с работами других учёных, занимавшихся колебаниями гравитационного поля Земли, периодичностями, связанными с Луной и Солнца, климатическими изменениями.

Поскольку временной ряд для J₂ ранее не исследовался с помощью сингулярного спектрального анализа, в нашем исследовании удалось подтвердить ранее выявленные закономерности и обнаружить новые. Важным результатом является обнаружение тренда, с возможным периодом около 60 лет, который также может быть связан с динамикой Чандлеровского колебания полюса Земли и изменениями длительности земных суток. Так как эти процессы имеют период в 60 лет, то можно предполагать, что между ними имеется взаимосвязь.

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

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

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

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

Список источников

1. Жарков В.Н. Внутреннее строение Земли и планет. Элементарное введение в планетную и спутниковую геофизику. М.: Наука и образование, 2013.

2. Minkang Cheng and Byron D. Tapley A 33 Year Time History of the J2 Changes from SLR // ILRS Workshop, 2008

3. G.A. Krasinsky. Variations of the coefficient J₂ of geopotential, and the dynamical Love number from the analysis of laser ranging to LAGEOS 1 and LAGEOS 2. // Cornell University Library, 2011

4. Кашкин В.Б., Рублева Т.В. Применение сингулярного спектрального анализа для выделения слабо выраженных трендов // Известия ТПУ. 2007. №5.

5. Зотов Л.В., Бизуар К., Шам С.К. О возможной взаимосвязи вращения земли и изменений климата в последние 150 лет. // Труды Всероссийской астрометрической конференции «ПУЛКОВО - 2015»

6. NASA. Planetary Fact Sheet. URL: http://nssdc.gsfc.nasa.gov/planetary/factsheet/

(дата обращения 27.05.2016).

7. Beauvalet, L., V. Lainey, J.-E. Arlot, and R.P. Binzel. «Dynamical Parameter Determinations in Pluto's System.» Astronomy & Astrophysics 540 (April 2012): A65.

8. A. Aitta. Internal structure of Pluto and Charon with an iron core. // Cornell University Library, 2015

9. Кузнецов В.В. Физика Земли: Учебник-монография. Новосибирск, 2011

10. Long-term 30-day estimates of C20 from up to 8 SLR satellites. URL: ftp://ftp.csr.utexas.edu/pub/slr/degree_2/Long_term/C20_1976_2011.txt (дата обращения 27.05.2016).

11. Description for UT/CSR monthly C20 RL-05 time series from SLR. URL: ftp://ftp.csr.utexas.edu/pub/slr/degree_2/C20_RL05.txt

(дата обращения 27.05.2016).

12. Xiaoping W., Jim Ray b, Tonie van Dam. Geocenter motion and its geodetic and geophysical implications // Journal of Geodynamics, 2012.

13. L. Zotov, Ch. Bizouard, C.K. Shum. About possible interrelation between Earth rotation and Climate variability on a decadal time-scale? // Journal of Geodesy and Geodynamics, China, 2016

14. Last decade's slow-down in global warming enhanced by an unusual climate anomaly // JOINT Research center, 2014 URL: https://ec.europa.eu/jrc/en/news/climate-anomaly-causes-global-warming-slow-down (дата обращения 27.05.2016)

15. Keeling C.D. Whorf T.P. Decadal oscillations in global temperature and atmospheric carbon dioxide. // Climate Research Committee, National Research Council. Washington, DC. 1996.

16. Пантелеев В.Л. Физика Земли и планет: Курс лекций. М: Московский государственный университет им. М.В. Ломоносова, физический факультет, 2001.

17. Мюррей К., Дермотт С. Динамика Солнечной системы. М: Физматлит, 2010.

18. Холшевников К.В., Титов В.Б. Задача двух тел: Учеб. пособие. СПб: СПБГУ, 2007.

19. Зотов Л.В. Теория фильтрации и обработка временных рядов: Курс лекций. М: Физический факультет МГУ, 2010

20. Satellite Laser Ranging and Earth Science // NASA International Laser Ranging Service, 2009.

21. Миллер Н.О., Воротков М.В. Моделирование чандлеровского движения полюса. // Труды Всероссийской астрометрической конференции «ПУЛКОВО - 2015»

22. Сидоренков Н.С. О важности наблюдений за эффектами приливных колебаний скорости вращения Земли. URL: http://meteoweb.ru/articles/sidorenkov.pdf (дата обращения 27.05.2016)

23. D.G. Kiryan, G.V. Kiryan. Moon's perigee mass as a missing component of the Earth's precession-nutation theory. // Cornell University Library, 2013

24. Зотов Л., Фролова Н., Шам С.К., Гравитационные аномалии в бассейнах крупных рек России, Природа, РАН, N5 2016, стр. 3

Размещено на Allbest.ru