Космические исследования, 2022, T. 60, № 2, стр. 134-150

ВВЕДЕНИЕ

Данная статья является продолжением работ [1, 2]. В [1] проведена повторная обработка магнитных измерений, выполненных в 1999 г. на спутнике Фотон-12. Первоначальная обработка была сделана вскоре после полета [3] с использованием детальной модели вращательного движения этого спутника. Модель, принятая в [1], проще модели [3], но сравнима с ней по точности при реконструкции движения с не очень малой угловой скоростью. Модель [1] позволила реконструировать эволюцию фактического движения Фотон-а-12 с возрастающей угловой скоростью. Ниже используется еще более простая модель, в которой реальная орбита спутника (высота апогея 380 км, высота перигея 220 км), заменена круговой орбитой. Эта замена позволила описать вращательное движение спутника автономными дифференциальными уравнениями, которые в случае установившегося движения с достаточно большой угловой скоростью (результат указанной выше эволюции) можно свести к обобщенно-консервативной системе. Такое сведение позволило приближенно, но достаточно точно описать реальное движение спутника на отрезках времени до 3.5 ч периодическими решениями, продолженными из периодических решений Ляпунова. Подходящие для этой цели периодические решения были изучены в [2]. Там же отмечалось их сходство с фактически вращательным движением Фотон-12. Однако такое сходство не было исследовано детально. Ниже такое исследование проводится.

УРАВНЕНИЯ ВРАЩАТЕЛЬНОГО ДВИЖЕНИЯ СПУТНИКА

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

Система $O{{x}_{1}}{{x}_{2}}{{x}_{3}}$ образована главными центральными осями инерции спутника. Точка $O$ – центр масс спутника, ось $O{{x}_{1}}$ – ось материальной симметрии спутника. Эта ось близка к продольной оси спутника и направлена от спускаемого аппарата к приборному отсеку. Момент инерции спутника относительно оси $O{{x}_{1}}$ обозначим ${{I}_{1}},$ равные моменты инерции относительно осей $O{{x}_{2}}$ и $O{{x}_{3}}$ обозначим ${{I}_{2}}.$

Вспомогательная система координат $O{{y}_{1}}{{y}_{2}}{{y}_{3}}$ служит для записи уравнений вращательного движения спутника. Ось $O{{y}_{1}}$ совпадает с осью $O{{x}_{1}};$ оси $O{{x}_{2}},$ $O{{x}_{3}}$ получаются из осей $O{{y}_{2}},$ $O{{y}_{3}}$ поворотом системы $O{{y}_{1}}{{y}_{2}}{{y}_{3}}$ на угол $\chi $ вокруг оси $O{{y}_{1}}.$ Кинематическая связь между системами $O{{x}_{1}}{{x}_{2}}{{x}_{3}}$ и $O{{y}_{1}}{{y}_{2}}{{y}_{3}}$ задается условием, что проекция абсолютной угловой скорости второй из них на ось $O{{y}_{1}}$ равна нулю. Проекции этой угловой скорости на оси $O{{y}_{2}},$ $O{{y}_{3}}$ обозначим ${{w}_{2}},$ ${{w}_{3}}.$ Пусть абсолютная угловая скорость спутника ${\mathbf{\omega }}$ имеет в системе $O{{x}_{1}}{{x}_{2}}{{x}_{3}}$ компоненты $({{\omega }_{1}},{{\omega }_{2}},{{\omega }_{3}}).$ Тогда $\dot {\chi } = {{\omega }_{1}}$ и

Здесь и ниже точкой обозначается дифференцирование по времени t.

В приборной системе координат $O{{z}_{1}}{{z}_{2}}{{z}_{3}}$ интерпретируются данные измерений бортовых магнитометров. Эту систему можно перевести в систему $O{{x}_{1}}{{x}_{2}}{{x}_{3}}$ двумя последовательными поворотами. Первый поворот выполняется на угол ${{\alpha }_{c}}$ вокруг оси $O{{z}_{2}},$ второй поворот – на угол ${{\beta }_{c}}$ – выполняется вокруг оси $O{{z}_{3}},$ получившийся после первого поворота. В общем случае, чтобы задать положение одной системы координат относительно другой, необходимы три угла. В данном случае можно было бы ввести еще угол поворота приборной системы вокруг ее оси $O{{z}_{1}},$ получившейся после первых двух поворотов. Однако поскольку направление одной из осей $O{{x}_{2}},$ $O{{x}_{3}}$ можно выбирать произвольно, третий угол удобно принять равным нулю, фиксировав тем самым положение системы $O{{x}_{1}}{{x}_{2}}{{x}_{3}}$ относительно системы $O{{z}_{1}}{{z}_{2}}{{z}_{3}}.$ Матрицу перехода между этими системами координат обозначим $\left\| {{{b}_{{ij}}}} \right\|_{{i,j = 1}}^{3},$ где ${{b}_{{ij}}}$ – косинус угла между осями $O{{z}_{i}}$ и $O{{x}_{j}}.$ Элементы этой матрицы выражаются через углы ${{\alpha }_{c}}$ и ${{\beta }_{c}}$ по формулам

Вращательное движение спутника изучается относительно орбитальной системы координат $O{{X}_{1}}{{X}_{2}}{{X}_{3}}.$ Ее оси $O{{X}_{1}}$ и $O{{X}_{3}}$ направлены по геоцентрическим скорости и радиус-вектору точки $O$. Матрицу перехода от системы $O{{y}_{1}}{{y}_{2}}{{y}_{3}}$ к орбитальой системе обозначим $\left\| {{{a}_{{ij}}}} \right\|_{{i,j = 1}}^{3},$ ${{a}_{{ij}}}$ – косинус угла между осями $O{{X}_{i}}$ и $O{{y}_{j}}.$ Элементы этой матрицы будем выражать в функции углов $\psi ,$ $\theta $ и $\delta ,$ которые введем так, чтобы система $O{{X}_{1}}{{X}_{2}}{{X}_{3}}$ переводилась в систему $O{{y}_{1}}{{y}_{2}}{{y}_{3}}$ тремя последовательными поворотами: 1) на угол $\psi $ вокруг оси $O{{X}_{3}},$ 2) на угол $\theta $ вокруг новой оси $O{{X}_{2}},$ 3) на угол $\delta $ вокруг новой оси $O{{X}_{1}},$ совпадающей с осью $O{{y}_{1}}.$ Таким образом, $\theta $ – угол между осью $O{{y}_{1}}$ и плоскостью $O{{X}_{1}}{{X}_{2}},$ $\psi $ – угол между проекцией оси $O{{y}_{1}}$ на плоскость $O{{X}_{1}}{{X}_{2}}$ и осью $O{{X}_{1}}.$ Эти два угла задают направление оси $O{{y}_{1}}$ в орбитальной системе координат. Имеют место соотношения

Система уравнений вращательного движения спутника образована динамическими уравнениями Эйлера для угловых скоростей ${{w}_{2}},$ ${{w}_{3}}$ и кинематическими уравнениями Пуассона для первой и третьей строк матриц $\left\| {{{a}_{{ij}}}} \right\|.$ В уравнениях Эйлера учитываются гравитационный и восстанавливающий аэродинамический моменты, а также постоянный момент вдоль оси $O{{x}_{1}}.$ При вычислении аэродинамического момента атмосфера считается неподвижной в абсолютном пространстве, ее плотность вдоль орбиты постоянна, внешняя оболочка спутника принимается сферой с центром на оси $O{{x}_{1}}.$ При сделанных допущениях аэродинамический момент характеризуется одним скалярным параметром. Система уравнений вращательного движения имеет вид [4, 5]

Здесь ${{\omega }_{0}}$ – орбитальная частота, $p$ – аэродинамический параметр, $\Omega $ и $\varepsilon $ – постоянные величины. В (2) использован явный вид решения одного из уравнений Эйлера ${{\dot {\omega }}_{1}} = \varepsilon $ с начальным условием ${{\omega }_{1}}({{t}_{0}}) = \Omega .$ Выбор ${{t}_{0}}$ будет указан ниже.

При численном интегрировании уравнений (2) единицей измерения времени служит 10³ c, единицы измерения других величин: $[{{\omega }_{i}}] = [{{w}_{i}}] = {{10}^{{ - 3}}}\,\,{{{\text{c}}}^{{ - 1}}},$ $[\varepsilon ] = [p] = {{10}^{{ - 6}}}\,\,{{{\text{c}}}^{{ - 2}}}.$ Переменные ${{a}_{{1i}}}$ и ${{a}_{{3i}}}$ зависимы, они связаны условиями ортогональности матрицы $\left\| {{{a}_{{ij}}}} \right\|.$ По этой причине начальные условия ${{a}_{{1i}}}({{t}_{0}})$ и ${{a}_{{3i}}}({{t}_{0}})$ выражаются через углы $\psi ,$ $\theta $ и $\delta .$ Элементы второй строки матрицы $\left\| {{{a}_{{ij}}}} \right\|$ вычисляются как векторное произведение ее третьей и первой строк. Формулы (1) и соотношение $\chi = \Omega (t - {{t}_{0}}) + {{\varepsilon {{{(t - {{t}_{0}})}}^{2}}} \mathord{\left/ {\vphantom {{\varepsilon {{{(t - {{t}_{0}})}}^{2}}} 2}} \right. \kern-0em} 2}$ позволяют найти функции ${{\omega }_{2}}(t),$ ${{\omega }_{3}}(t)$ и движение системы $O{{x}_{1}}{{x}_{2}}{{x}_{3}},$ решая уравнения (2).

Параметр $\lambda $ известен: $\lambda \approx 0.24.$ Тем не менее, он и параметры $p$ и $\varepsilon $ определяются из обработки данных измерений наряду с неизвестными начальными условиями движения спутника, т.е. служат параметрами согласования.

Уравнения (2) получены из аналогичных уравнений вращательного движения [1] переходом к круговой орбите и использованием орбитальной системы координат в качестве системы, относительно которой рассматривается движение спутника. Для компенсации всех сделанных упрощений реконструировались движения спутника, в которых компонента угловой скорости ${{\omega }_{1}}$ достаточно велика. Обработка измерений с помощью уравнений (2) выполнялась по схеме [1] и дала близкие результаты.

РЕКОНСТРУКЦИЯ ВРАЩАТЕЛЬНОГО ДВИЖЕНИЯ СПУТНИКА ПО МАГНИТНЫМ ИЗМЕРЕНИЯМ

На борту Фотона-12 находилась аппаратура “Мираж” с несколькими трехкомпонентными магнитометрами. Поскольку движение спутника было неуправляемым, данные измерений этой аппаратуры и уравнения (2) можно использовать для определения фактического вращательного движения спутника по обычным статистическим методикам. Методика, использованная ниже, состоит в следующем [1, 4]. По измерениям, выполненным на некотором отрезке времени ${{t}_{0}} \leqslant t \leqslant {{t}_{0}} + T,$ строились функции ${{\hat {h}}_{i}}(t),$ которые задавали на этом отрезке компоненты вектора местной напряженности магнитного поля в системе координат $O{{z}_{1}}{{z}_{2}}{{z}_{3}}.$ Среднеквадратичные ошибки аппроксимации не превышали 200$\gamma $ ($1\gamma = {{10}^{{ - 5}}}$ Э). Затем вычислялись псевдоизмерения $h_{i}^{{(n)}} = {{\hat {h}}_{i}}({{t}_{n}}),$ ${{t}_{n}} = {{t}_{0}} + {{nT} \mathord{\left/ {\vphantom {{nT} N}} \right. \kern-0em} N},$ где $n = 0,1,2,...N.$ Обычно было $T = 100{\kern 1pt} - {\kern 1pt} 300$ мин, ${T \mathord{\left/ {\vphantom {T N}} \right. \kern-0em} N} \approx 1$ мин. Псевдоизмерения служили исходной информацией для отыскания решений уравнений (2), описывающих фактическое движение спутника.

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

Здесь ${{\Delta }_{i}}$ – оценки постоянных смещений в псевдоизмерениях, $h_{j}^{ \circ }(t)$ – компоненты напряженности МПЗ в точке $O$ в системе координат $O{{x}_{1}}{{x}_{2}}{{x}_{3}},$ рассчитываемые с помощью модели IGRF2005.

Чтобы использовать уравнения (2) для вычисления функционала (3), необходимо задать связь системы $O{{X}_{1}}{{X}_{2}}{{X}_{3}}$ с гринвичской системой координат, в которой задана реальная орбита спутника и вычисляется напряженность магнитного поля Земли. С этой целью фактическая орбита спутника на отрезке ${{t}_{0}} \leqslant t \leqslant {{t}_{0}} + T$ аппроксимировалась круговой орбитой. Аппроксимация строилась методом наименьших квадратом по достаточно точным значениям реального фазового вектора центра масс спутника, заданным на равномерной сетке с шагом 3 мин. Круговая орбита задавалась 5-ю элементами: радиусом, средним движением (орбитальной частотой ${{\omega }_{0}}$), начальным значением аргумента широты, долготой восходящего узла и наклонением. Вдоль круговой орбиты вычислялись компоненты напряженности МПЗ в системе координат $O{{X}_{1}}{{X}_{2}}{{X}_{3}}$ на моменты ${{t}_{n}}.$ Эти компоненты использовались при многократном вычислении функционала (3) в процессе его минимизации. С помощью решения уравнений (2) они пересчитывались в систему $O{{y}_{1}}{{y}_{2}}{{y}_{3}}$ и затем в систему $O{{x}_{1}}{{x}_{2}}{{x}_{3}}.$ В качестве ${{t}_{0}}$ в уравнениях (2) всегда использовалась начальная точка обрабатываемого отрезка данных.

Оценку приемлемости сделанного перехода к круговой орбите, а также оценку ожидаемой точности аппроксимации псевдоизмерений экстремалью функционала (3) дает минимизация функции

по $\kappa $ и $\Delta _{i}^{'}.$ Здесь $\kappa \approx 1$ – масштабирующий коэффициент, $\Delta _{i}^{'}$ – оценки постоянных смещений в псевдоизмерениях. Первый радикал в формуле для $\Psi $ – скорректированный модуль измеренной напряженности магнитного поля на спутнике в момент ${{t}_{n}},$ второй радикал – расчетный модуль напряженности МПЗ в точке $O$ в тот же момент. Функция $\Psi $ характеризует близость модулей измеренной и расчетной напряженности магнитного поля на отрезке ${{t}_{0}} \leqslant t \leqslant {{t}_{0}} + T.$ Для расчета $\Psi $ надо знать только орбитальное движение спутника. Минимальное значение ${{\Psi }_{{\min }}}$ этой функции дает предварительную оценку точности магнитных измерений.

При $T = 210$ мин для большинства отрезков псевдоизмерений среднеквадратичное отклонение ${{\sigma }_{ * }} = \sqrt {{{{{\Psi }_{{\min }}}} \mathord{\left/ {\vphantom {{{{\Psi }_{{\min }}}} {(N - 3)}}} \right. \kern-0em} {(N - 3)}}} $ лежит в пределах 1600–2000γ. При расчете ${{\Psi }_{{\min }}}$ с использованием реальной орбиты ${{\sigma }_{ * }}$ меньше и лежит в пределах 1400–1600γ. Два примера минимизации функции $\Psi $ приведены на рис. 1. В верхней части рисунка сравниваются модули вектора напряженности МПЗ, рассчитанные по скорректированным псевдоизмерениям и по модели IGRF. Расчет по псевдоизмерениям представлен маркерами, расчет по модели IGRF – сплошными линиями. Графики в нижней части рисунка характеризуют отклонения псевдоизмерений от модели. Эти графики – ломаные, вершины которых отвечают ошибкам аппроксимации. В подписи к рисунку указаны значения ${{\sigma }_{ * }}$; рядом с ними в скобках указаны значения ${{\sigma }_{ * }}$, полученные с использованием реального орбитального движения. В функционал (3) подставлялись псевдоизмерения, масштаб которых был скорректирован множителем $\kappa .$ Значения этого множителя уточнялись для каждого обрабатываемого отрезка данных.

Функционал (3) минимизировался по 11 величинам: начальным условиям решения системы (2) $\psi ({{t}_{0}}),$ $\theta ({{t}_{0}}),$ $\delta ({{t}_{0}}),$ $\Omega ,$ ${{w}_{2}}({{t}_{0}}),$ ${{w}_{3}}({{t}_{0}})$ и параметрам $\lambda ,$ $p,$ $\varepsilon ,$ ${{\alpha }_{c}},$ ${{\beta }_{c}}.$ Заключительный этап минимизации выполнялся методом Гаусса–Ньютона. Точность аппроксимации псевдоизмерений и разброс в оцениваемых величинах характеризовались соответствующими стандартными отклонениями. Стандартные отклонения рассчитывались в предположении, что ошибки в псевдоизмерениях $h_{i}^{{(n)}}$ некоррелированы и имеют одинаковые дисперсии, средние значения ошибок в псевдоизмерениях с одинаковым нижним индексом $i$ одинаковы (величины ${{\Delta }_{i}}$ в (3) – оценки этих средних значений). Стандартные отклонения вычислялись так. Пусть ${{\Phi }_{{\min }}}$ – значение функционала (3) в точке минимума, $C$ – матрица системы нормальных уравнений метода Гаусса–Ньютона в этой точке (матрица $2C$ приблизительно равна матрице квадратичной формы ${{d}^{2}}\Phi $ в точке минимума $\Phi $). Тогда дисперсия ошибок в псевдоизмерениях оценивается величиной $\sigma _{H}^{2} = {{{{\Phi }_{{\min }}}} \mathord{\left/ {\vphantom {{{{\Phi }_{{\min }}}} {(3N - 11)}}} \right. \kern-0em} {(3N - 11)}}.$ Стандартные отклонения оцениваемых величин равны квадратным корням из соответствующих диагональных элементов матрицы $\sigma _{H}^{2}{{C}^{{ - 1}}}.$ Эти стандартные обозначим ${{\sigma }_{\psi }},$ ${{\sigma }_{\theta }},$ ${{\sigma }_{\delta }},$ ${{\sigma }_{\Omega }},$ ${{\sigma }_{{w2}}},$ ${{\sigma }_{{w3}}},$ ${{\sigma }_{\lambda }},$ ${{\sigma }_{p}},$ ${{\sigma }_{\varepsilon }},$ ${{\sigma }_{{\alpha c}}},$ ${{\sigma }_{{\beta c}}}.$

Движение спутника Фотон-12 реконструировано на 25 временных интервалах [5]. Некоторые результаты реконструкции, относящиеся к интервалам 10–25 приведены в табл. 1, реконструкции движения на интервалах 14 и 17 представлены на рис. 2–4. В подписях к рисункам и в таблице использовано Всемирное координированное время (UTC). В табл. 1 приведены начальные точки интервалов ${{t}_{0}}$ и их фактическая длина T. Номинальная длина $T = 210$ мин выбрана таким образом, чтобы получаемые значения ${{\sigma }_{H}}$ не более чем в 2 раза превышали указанные выше типичные значения ${{\sigma }_{ * }}.$ Иными словами, чтобы принятая модель движения была достаточно адекватна. Поскольку принятая модель является упрощенной, с ее помощью реконструировались движения спутника с заметной угловой скоростью, которые возникли через несколько (~2) суток после начала неуправляемого движения. Табл. 1 содержит данные по интервалам, на которых угловая скорость достаточно большая.

Рис. 2 иллюстрирует типичное качество аппроксимации псевдоизмерений. Здесь сплошными линиями изображены графики функций ${{h}_{i}}(t)$ при ${{t}_{0}} \leqslant t \leqslant {{t}_{0}} + T,$ маркеры указывают точки $\left( {{{t}_{n}},h_{i}^{{(n)}} - {{\Delta }_{i}}} \right),$ $n = 0,1,2,...,N.$ Количественно аппроксимация псевдоизмерений характеризуется стандартными отклонениями ${{\sigma }_{H}},$ значения которых приведены в табл. 1 и в подписи к рисунку. Значения ${{\sigma }_{H}}$ в табл. 1 несколько выше, чем в [1]. Тем не менее, достигнутая точность достаточна для целей данной работы.

На рис. 3 помещены графики компонент угловой скорости ${{w}_{2}}$ и ${{w}_{3}}.$ Компонента угловой скорости ${{\omega }_{1}}$ в движении на рис. 3а монотонно убывает, а в движении на рис. 3б монотонно возрастает. Пределы изменения этой переменной указаны в подписи к рисунку.

По мере увеличения компоненты ${{\omega }_{1}}$ движение спутника становилось все больше похоже на регулярную прецессию Эйлера осесимметричного твердого тела. Формирование регулярной прецессии с медленно возрастающей угловой скоростью ${{\omega }_{1}}$ и мало меняющимся значением ${{\omega }_{ \bot }} = \sqrt {w_{2}^{2} + w_{3}^{2}} $ завершилось после нескольких суток полета.

Точная регулярная прецессия Эйлера может иметь место лишь в случае, когда спутник осесимметричен и главный момент приложенных к нему внешних сил равен нулю. Тогда величины ${{\omega }_{1}}$ и $\omega {}_{ \bot }$ остаются неизменными во время движения. Для используемой модели первое условие выполнено, а второе – нет. По этой причине можно говорить лишь о движениях, близких к регулярной прецессии. Такие движения удобно характеризовать величинами

и определяемыми аналогичными формулами величинами ${{\bar {\omega }}_{1}},$ $\delta {{\omega }_{1}}.$ Для последних имеют место формулы Среднеквадратичные отклонения $\delta {{\omega }_{1}}$ и $\delta {{\omega }_{ \bot }}$ характеризуют близость движения спутника к регулярной прецессии с параметрами ${{\bar {\omega }}_{1}}$ и ${{\bar {\omega }}_{ \bot }}.$ Значения величин ${{\bar {\omega }}_{1}},$ $\delta {{\omega }_{1}},$ ${{\bar {\omega }}_{ \bot }}$ и $\delta {{\omega }_{ \bot }}$ указаны в табл. 1.

На рис. 4 приведены графики зависимости от времени углов $\psi $ и $\theta ,$ задающих положение оси $O{{x}_{1}}$ относительно координат $O{{X}_{1}}{{X}_{2}}{{X}_{3}}.$ Эти углы рассчитывались с использованием формул

Движение спутника по углам установилось только к интервалу 7 – ось $O{{x}_{1}}$ перестала пересекать плоскость орбиты и весь дальнейшей полет оставалась в полупространстве ${{Z}_{2}} > 0.$ Движение на рис. 3, 4 уже явно похоже на регулярную прецессию Эйлера.

Дополнительную информацию о движении спутника дает рис. 5. Здесь для двух рассмотренных выше реконструкций приведены графики угла $\alpha $ между осью $O{{x}_{1}}$ и вектором кинетического момента спутника в его движении относительно центра масс и угла $\rho $ между этим вектором и осью $O{{X}_{2}}.$ В регулярной прецессии Эйлера $\alpha = {\text{const}}{\text{.}}$ В примерах на рис. 5 этот угол меняется в достаточно узких пределах вблизи значений ${{\alpha }_{ * }} = {\text{arctg(}}{{{{{\bar {\omega }}}_{ \bot }}} \mathord{\left/ {\vphantom {{{{{\bar {\omega }}}_{ \bot }}} {\lambda {{{\bar {\omega }}}_{1}}}}} \right. \kern-0em} {\lambda {{{\bar {\omega }}}_{1}}}}{\text{),}}$ указанных в подписи к рисунку. Вектор кинетического момента колеблется вблизи оси $O{{X}_{2}}$ (см. графики угла $\rho $), угловая скорость которой менее 5 град/сут. Движения такого типа имели место 15–19.IX, в частности, на интервалах 10 и 14–18. Оценки параметров $p,$ $\lambda ,$ $\varepsilon ,$ ${{\alpha }_{c}},$ ${{\beta }_{c}}$ и их стандартные отклонения для этих интервалов приведены в табл. 2. Здесь углы выражены в радианах, величины $p$ и $\varepsilon $ – в единицах ${{10}^{{ - 6}}}\,{{{\text{c}}}^{{ - 2}}}.$ Данные этой таблицы характерны для всех интервалов 1–25 [5]. Стандартные отклонения начальных условий реконструкций на интервалах 1–25, выраженные в радианах и 0.001 с^–1, лежат в диапазонах [5]: ${{\sigma }_{\delta }} = 0.01{\kern 1pt} - {\kern 1pt} 0.025,$ ${{\sigma }_{{\theta ,\psi }}} = 0.01{\kern 1pt} - {\kern 1pt} 0.02,$ ${{\sigma }_{\Omega }} = 0.01{\kern 1pt} - {\kern 1pt} 0.025,$ ${{\sigma }_{{w2,w3}}} = 0.02{\kern 1pt} - {\kern 1pt} 0.08.$ Приведенные примеры демонстрируют менее точную, чем в [1, 3], аппроксимацию псевдоизмерений. Значения ${{\sigma }_{H}}$ в табл. 1 несколько выше, чем в [1, 3]. Тем не менее, реконструкции данной работы передают все особенности вращательного движения спутника.

Интервалы 10 и 14–18 представляют особый интерес для дальнейшего, поскольку движения на них допускают приемлемую аппроксимацию периодическими движениями обобщенно-консервативной механической системы.

ОБОБЩЕННО-КОНСЕРВАТИВНАЯ МОДЕЛЬ ВРАЩАТЕЛЬНОГО ДВИЖЕНИЯ СПУТНИКА

Как видно из табл. 1, значения ${{\bar {\omega }}_{1}}$ со временем стабилизировались. По-видимому, возрастал не учитываемый в уравнениях (2) момент сопротивления, пропорциональный ${{\omega }_{1}}$ [1]. На интервалах из той части табл. 1, где значения $\varepsilon $ малы, для реконструкции движения спутника можно использовать решения уравнений (2) при $\varepsilon = 0.$ Графики, иллюстрирующие такую реконструкцию на интервалах 10 и 14–18, приведены в [5]. Они очень похожи на графики, полученные с помощью исходной системы (2), но имеют увеличенные значения ${{\sigma }_{H}}.$ Эти значения на указанных интервалах приведены в табл. 3 вместе с оценками величин $\Omega ,$ $p,$ $\lambda ,$ ${{\alpha }_{c}},$ ${{\beta }_{c}}$ и их стандартными отклонениями. В этой таблице $[\Omega ] = 0.001\,\,{{{\text{c}}}^{{ - 1}}},$ $[p] = {{10}^{{ - 6}}}\,\,{{{\text{c}}}^{{ - 2}}}.$ Сравнение данных в табл. 2 и 3 позволяет понять, что дало сделанное огрубление модели. Найденные при новой обработке значения $\Omega $ весьма близки соответствующим значениям ${{\bar {\omega }}_{1}}$ из табл. 1.

Ниже в этом разделе полагаем в (2) $\varepsilon = 0,$ ${{\omega }_{1}} = \Omega $ и считаем $\Omega $ параметром. Уравнения (2) дополним уравнениями Пуассона для компонент орта оси $O{{X}_{2}}$ в системе $O{{x}_{1}}{{x}_{2}}{{x}_{3}}{\text{:}}$

и получившуюся систему обозначим (2a). Уравнения (2a) допускают обобщенный интеграл энергии т.е. спутник представляет собой обобщенно-консервативную механическую систему. Уравнения (2a) имеют также семейства частных решений, в которых ${{a}_{{i1}}}$ – постоянные величины, связанные соотношениями а остальные элементы матрицы $\left\| {{{a}_{{ij}}}} \right\|$ определяются соотношениями Пусть соотношения (4) выполнены. Уравнения (5) для переменных ${{a}_{{22}}}$ и ${{a}_{{23}}}$ имеют решение ${{a}_{{22}}} = \sqrt {1 - a_{{21}}^{2}} \cos [{{a}_{{21}}}{{\omega }_{0}}(t - {{t}_{0}})],$ ${{a}_{{23}}} = \sqrt {1 - a_{{21}}^{2}} \times $ $ \times \,\sin [{{a}_{{21}}}{{\omega }_{0}}(t - {{t}_{0}})],$ где ${{t}_{0}}$ – произвольная постоянная; переменные ${{a}_{{12}}},$ ${{a}_{{13}}},$ ${{a}_{{32}}},$ ${{a}_{{33}}}$ находятся по формулам (5). Орт оси Ox₁ в системе OX₁X₂X₃ имеет компоненты ${{a}_{{i1}}},$ поэтому решения (4), (5) описывают покой оси Ox₁ в этой системе. При $p = 0$ эти решения совпадают с известными решениями, называемыми конической, цилиндрической и гиперболоидальной прецессиями [6]. Их сложный вид и неожиданный период связаны со способом введения системы координат Oy₁y₂y₃. Эта система, удобная в задаче обработки магнитных измерений, плохо подходит для описания движения спутника относительно орбитальной системы координат. Оказалось, что на интервалах 10, 14–18 для реконструкции можно использовать периодические движения специального вида. Но чтобы прийти к ним, уравнения (2a) надо заменить эквивалентными уравнениями в других переменных.

Для описания движений спутника, в которых ось Ox₁ отклоняется от оси $O{{X}_{2}}$ менее чем на 90°, положение системы Ox₁x₂x₃ относительно системы OX₁X₂X₃ удобно задавать введенными выше углами $\psi ,$ $\theta $ и третьим углом $\varphi ,$ поворот на который вокруг оси $O{{x}_{1}}$ завершает преобразование системы OX₁X₂X₃ в систему Ox₁x₂x₃. Преобразование OX₁X₂X₃ → Ox₁x₂x₃ можно представить в виде суперпозиции преобразований OX₁X₂X₃ → Oy₁y₂y₃ и Oy₁y₂y₃ → Ox₁x₂x₃. Первое преобразование задается углами $\psi ,$ $\theta $ и $\delta ,$ второе преобразование – углом $\chi .$ Повороты на углы $\delta $ и $\chi $ выполняются вокруг оси $O{{x}_{1}} = O{{y}_{1}}$ в одном направлении, поэтому $\varphi = \delta + \chi .$

Если интересоваться только движением оси симметрии спутника $O{{x}_{1}},$ то удобно использовать углы $\psi ,$ $\theta $ и комбинации угловых скоростей

где Такие переменные были использованы в [2]. Примеры графиков угловых скоростей $\dot {\delta },$ ${{\Omega }_{2}},$ ${{\Omega }_{3}}$ в реконструкциях на интервалах 14 и 17 приведены на рис. 6. Величины ${{\Omega }_{2}},$ ${{\Omega }_{3}}$ – проекции угловой скорости спутника на оси Резаля, отвечающие углам $\psi $ и $\theta .$

ПЕРИОДИЧЕСКИЕ ДВИЖЕНИЯ ОСИ СИММЕТРИИ СПУТНИКА

Для переменных $\psi ,$ $\theta ,$ ${{\Omega }_{2}},$ ${{\Omega }_{3}}$ справедливы уравнения [2]

Выраженные через $\psi ,$ $\theta ,$ ${{\Omega }_{2}}$ и ${{\Omega }_{3}}$ решения системы (2a) на интервалах 10 и 14–18 похожи на периодические решения системы (6), изучавшиеся в [2]. Решения на этих интервалах отличаются от решений на остальных интервалах тем, что кинетический момент движения спутника относительно центра масс мало отклоняется от нормали к плоскости орбиты.

Система (6) инвариантна относительно преобразования переменных $t \to - t,$ $\theta \to - \theta ,$ ${{\Omega }_{3}} \to - {{\Omega }_{3}},$ поэтому можно искать ее решения, в которых переменные $\theta $ и ${{\Omega }_{3}}$ – нечетные, а переменные $\psi $ и ${{\Omega }_{2}}$ – четные функции времени. Периодические решения с такими свойствами называются симметричными. Симметричные периодические решения с периодом $T$ определяются краевыми условиями

В [2] построены симметричные периодические решения системы (6), продолженные из решений Ляпунова, существующих в окрестности ее стационарного решения $\theta = 0,$ $\psi = \psi (p),$ $\left| p \right| \leqslant 0.5,$ $\psi (0) = \pi {\text{/}}2.$ При $p = 0$ это стационарное решение переходит в цилиндрическую прецессию. Существует два семейства таких периодических решений, называемых короткопериодическими и длиннопериодическими.

Для аппроксимации использовались короткопериодические решения. Параметры системы (6) $p,$ $\lambda ,$ $\Omega $ для каждого аппроксимируемого интервала определялись в результате обработки магнитных измерений (табл. 3). Значение периода $T$ находились в результате спектрального анализа решений системы (2a), аппроксимирующих магнитные измерения на исследуемом интервале. Переменные системы (2a) в этом решении указанным выше способом пересчитывались в функции θ(t), $\psi (t),$ ${{\Omega }_{2}}(t),$ ${{\Omega }_{3}}(t),$ и для каждой из этих функций определялся пробный период. Опишем определение пробного периода на примере функции θ(t). Эта функция (как и остальные из приведенного набора) вычислялась на равномерной временной сетке $\left\{ {t_{m}^{'}} \right\}_{{m = 0}}^{M},$ $t_{m}^{'} = {{t}_{0}} + mh,$ $h = 20$ с, $(M - 1)h < T \leqslant Mh.$ Строилась ее аппроксимация выражением

где ${{a}_{0}},$ $a,$ $b$ и $f$ – параметры. Значения параметров искались методом наименьших квадратов посредством минимизации выражения Минимизация выполнялась в два этапа. Сначала в результате решения ряда одинаковых линейных задач наименьших квадратов вычислялись значения функции в узлах достаточно мелкой равномерной сетки на отрезке $0 \leqslant f \leqslant {{10}^{{ - 3}}}$ Гц. Затем на этой сетке находилась точка ${{f}_{{\min }}} = \arg \min {{\Psi }_{1}}(f).$ Пробный период определялся по формуле $T = {1 \mathord{\left/ {\vphantom {1 {{{f}_{{\min }}}}}} \right. \kern-0em} {{{f}_{{\min }}}}}.$ Если функция ${{\Psi }_{1}}(f)\,$ имеет несколько значимых минимумов, то исследуемая переменная θ(t) содержит несколько составляющих вида (8). В этом случае периодическая аппроксимация может оказаться невозможной.

На интервалах 10 и 14–18 функции ${{\Psi }_{1}}(f)$ для переменных θ(t), $\psi (t),$ ${{\Omega }_{2}}(t),$ ${{\Omega }_{3}}(t)$ имели только один значимый минимум, причем значения $\min {{\Psi }_{1}}$ отличались от значений этих функций вне малой окрестности точки ${{f}_{{\min }}}$ (они близки $\max {{\Psi }_{1}}$) в 10 и более раз [5].

Пробные периоды, найденные по функциям θ(t), $\psi (t),$ ${{\Omega }_{2}}(t)$ и ${{\Omega }_{3}}(t),$ совпадали достаточно точно. Так как в движениях на интервалах 10 и 14–18 наиболее интенсивно изменялись переменные ${{\Omega }_{2}}$ и ${{\Omega }_{3}},$ пробный период определялся по периодограммам функций ${{\Omega }_{2}}(t)$ и ${{\Omega }_{3}}(t).$ Для них значения ${{f}_{{\min }}}$ попадали в одни и те же узлы сетки с шагом $2 \cdot {{10}^{{ - 7}}}$ Гц. Найденное значение пробного периода использовалось для построения аппроксимирующего решения краевой задачи (6), (7). Эта задача решалась методом пристрелки (см. [2]). Начальным приближением неизвестных начальных условий $\psi (0)$ и ${{\Omega }_{2}}(0)$ служили амплитуды соответствующих выражений (8) при $f = {{f}_{{\min }}}.$ Можно также использовать значения переменных $\psi $ и ${{\Omega }_{2}}$ в аппроксимируем решении в точке ${{t}_{ * }},$ определяемой соотношением $\theta ({{t}_{ * }}) \approx 0.$

Примеры построенных периодических решений приведены на рис. 7, 8 и в табл. 4. В таблице указаны параметры решений задачи (6), (7): $T{\text{/}}2,$ $\psi (0)$ и ${{\Omega }_{2}}(0)$, а также орбитальная частота ${{\omega }_{0}}$ на данном интервале и коэффициент $a$ в характеристическом уравнении ${{(\rho - 1)}^{2}}({{\rho }^{2}} - 2a\rho + 1) = 0,$ определяющем мультипликаторы периодического решения. Здесь единица измерения времени – 1000 с, единицы измерения угловой скорости – 0.001 c^–1. Как видим, все найденные периодические решения орбитально устойчивы в первом приближении. На рисунках приведены графики переменных θ(t), $\psi (t),$ ${{\Omega }_{2}}(t)$ и ${{\Omega }_{3}}(t)$ в периодических решениях на отрезке $0 \leqslant t \leqslant 2T$ и проекции орта оси $O{{x}_{1}}$ и орта кинетического момента спутника на плоскость $O{{X}_{1}}{{X}_{3}}.$ Проекция орта кинетического момента лежит внутри проекции орта оси $O{{x}_{1}}.$

Сопоставление движений, найденных обработкой магнитных измерений, и периодических решений выполнялось по следующей схеме. Переменные θ(t), $\psi (t),$ ${{\Omega }_{2}}(t)$ и ${{\Omega }_{3}}(t)$ в периодическом решении аппроксимировались дискретными рядами Фурье [7]. Четные функции разлагались в ряд по косинусам, нечетные – в ряд по синусам. В каждом аппроксимирующем выражении использовалось 200 гармоник. Построенные ряды обозначим, $\tilde {\theta }(t),$ $\tilde {\psi }(t),$ ${{\tilde {\Omega }}_{2}}(t)$ и ${{\tilde {\Omega }}_{3}}(t).$ Составлялось выражение

которое минимизировалось по $\tau $ на достаточно мелкой сетке. Найденные таким образом значения $\tau $ приведены в табл. 4. В левой части рис. 9–14 в единых системах координат приведены графики функций θ(t), $\psi (t),$ ${{\Omega }_{2}}(t)$ и ${{\Omega }_{3}}(t)$ и соответствующих аппроксимирующих выражений $\tilde {\theta }(t - {{t}_{0}} - \tau ),$ $\tilde {\psi }(t - {{t}_{0}} - \tau ),$ ${{\tilde {\Omega }}_{2}}(t - {{t}_{0}} - \tau ),$ ${{\tilde {\Omega }}_{3}}(t - {{t}_{0}} - \tau ).$ В правой части этих рисунков представлены графики разностей $\Delta \theta = \theta (t) - \tilde {\theta }(t - {{t}_{0}} - \tau ),$ и т.п. Все графики построены по значениям на сетке $\left\{ {t_{m}^{'}} \right\}_{{m = 0}}^{M}.$ Как видим, периодическая аппроксимация получилась достаточно точной. Судя по графикам, функции $\Delta \theta (t)$ и $\Delta \psi (t)$ содержат низкочастотную компоненту, которая, по-видимому, как-то связана с длиннопериодическими решениями Ляпунова [2].

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