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

Ниже описывается реализация метода определения массы станции, предложенного в [1]. Масса находится по данным измерений американского акселерометра MAMS, полученным во время коррекции обиты предварительно тарированными двигателями причаливания и ориентации (ДПО) транспортно грузового корабля (ТГК) «Прогресс». Использованы данные, полученные во время коррекции 20.08.2004. Те же данные и данные, полученные во время коррекции 26.08.2004, использованы для анализа низкочастотных колебаний корпуса станции, вызываемых работой ДПО. Определены собственные частоты двигателей. Описываемое исследование выполнено в рамках технологического эксперимента “Тензор”.

2. Методика определения массы станции по измерениям MAMS. Импульс, корректирующий орбиту станции, выдается вдоль оси  строительной системы координат Служебного модуля (СМ). Во время выдачи импульса эта ось направлена по трансверсали к орбите станции. Из теоремы об изменении количества движения станции в проекции на ось  следует соотношение

.

Здесь  – масса станции (вместе с пристыкованными к ней кораблями),  – время,  – кажущееся ускорение центра масс станции,  – число работающих двигателей,  – проекция на ось  тяги одного двигателя,  – компонента вдоль оси  главного вектора негравитационных внешних сил, действующих на станцию. Двигатели считаются одинаковыми, в  наибольший вклад вносят силы аэродинамического торможения.

Акселерометр MAMS позволяет измерить кажущееся ускорение  – одна из его чувствительных осей параллельна оси  строительной системы координат СМ. Пусть измеренные значения  известны для моментов времени  ,  – результаты измерений. Значения кусочно-постоянной функции  известны из телеметрии, силой  можно пренебречь. Тогда отношение  можно найти методом наименьших квадратов [2] из условия минимума выражения

,

,     ,     .

Оценка указанного отношения имеет вид

.                                                      (2.1)

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

,     .                             (2.2)

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

.                                                       (2.3)

Дифференциал этой формулы можно представить в виде

.

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

.                                                (2.4)

Здесь  – относительная ошибка в определении . Оценки величин  и  находились в результате тарировки тяги ДПО ТГК «Прогресс».

3. Тарировка тяги двигателей грузового корабля. Для тарировки тяги ДПО ТГК «Прогресс М – 42» на участке его автономного полета (участок дальнего сближения со станцией) на втором суточном витке 12.08.2004 был выдан тарировочный импульс. Импульс выдавался вдоль продольной оси корабля с помощью восьми двигателей ДПО в режиме без поддержания ориентации (контур ориентации был отключен), с выключением системы управления движением ТГК после выдачи импульса. Двигатели работали в непрерывном режиме.

На борту ТГК «Прогресс» размещается бортовой измеритель приращения кажущейся скорости (БИПС). Показания этого прибора представляют собой моменты времени, в которые приращение кажущейся скорости оказывается кратным 4 см/с. Для обработки показаний БИПС воспользуемся следующими уравнениями движения корабля во время выдачи импульса:

,     ,     ,     .

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

         Решение выписанных уравнений относительно  имеет вид

 .           (3.1)

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

На основании формулы (3.1) телеметрические данные БИПС сглаживались полиномами

                                  (3.2)

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

.

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

.                                                     (3.3)

Расчет  выполнен при  кг. Как видно из таблицы, стандартное отклонение  при переходе от  к  уменьшается более чем в пять раз, а при переходе от  к  менее чем в два раза. Минимальное значение  достигается при , второе по величине – при . Построенные аппроксимирующие полиномы в случае  не вполне согласуются с решением (3.1). Некоторые их старшие коэффициенты отрицательны, тогда как в разложении решения (3.1) все коэффициенты при степенях  положительны. Ниже при оценке массы станции используются значения , полученные при  и . Как видно из таблицы, это два крайних случая по значениям  и .

Оценим точность определения  по формуле (3.3). Дифференциал этой формулы имеет вид

.

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

.                                                 (3.4)

Здесь  – относительная ошибка в определении . Приняв , получим для  значения, указанные в табл. 1.

Оценка значения  была получена и от двигателистов. Длительный опыт эксплуатации кораблей “Союз” и “Прогресс” позволил получить обширный экспериментальный материал по зависимости тяги ДПО от давления в баках с топливом и от давления в системе наддува. Анализ телеметрической информации показал, что стабилизация давления в магистралях двигательной установки грузового корабля на участке выдачи корректирующего импульса (в составе станции) 20.08.2005 произошла на 160 секунде. До конца работы двигательной установки давления не практически не менялись. Таким образом, на данном интервале было  Н. Погрешность данного способа оценки тяги достаточно точно не известна. По оценкам листов , что связано в первую очередь с конструктивными особенностями системы наддува двигательной уста новки.

4. Оценка массы станции. Методика определения массы станции, описанная в п. 2, была реализована четырьмя способами. В первом способе в формулах (2.1), (2.2) использовались все данные измерений MAMS кажущегося ускорения , относящиеся к временному интервалу выдачи корректирующего импульса. График этих данных представлен на рис. 1а. График представляет собой ломаную, ординаты вершин которой суть данные измерений, а абсциссы образуют равномерную сетку с шагом  с. Во втором способе в формулы (2.1), (2.2) подставлялись сглаженные значения кажущегося ускорения. Сглаживание выполнялось так. Пусть   – данные измерений MAMS,  – натуральное число, . Сначала по этим данным стандартным методом наименьших квадратов строилось аппроксимирующее выражение

.

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

.

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

Третий и четвертый способы поясняются графиком функции  на рис. 1в. Напомним, эта кусочно-постоянная функция задает число ДПО, задействованных в момент  времени . Сравнение графиков на рис 1б и 1в показывает, что в формулах п. 2 целесообразно оставить лишь значения ускорения , для которых . Во-первых, суммарная тяга ДПО на этих участках наиболее стабильна; во-вторых, калибровка ДПО проводилась при одновременной работе именно восьми двигателей. В третьем способе в формулах п. 2 учитывались только низкочастотные данные MAMS, относящиеся к десяти достаточно продолжительным интервалам постоянства , на которых . В четвертом способе для тех же , что и в способе 3, использовались сглаженные значения .

Результаты определения массы станции по формулам (2.1), (2.3) для четырех описанных способов выбора величин  представлены в табл. 2, 3. Табл. 2 построена при  Н и , табл. 3 – при  Н и . Расчеты в обоих случаях выполнялись с одними и теми же величинами , поэтому значения  и  в этих таблицах одинаковы[1]. Указанная в таблицах относительная погрешность определения массы  рассчитана по формулам (2.4), (3.4). По результатам проведенной перед описываемым экспериментом инвентаризации проектное значение массы станции (вместе с ТГК «Прогресс») на момент начала выдачи корректирующего импульса 20.08. 2004 г. составляло  кг. Если считать это значение точным, то фактическую относительную погрешность определения массы станции можно оценить, как  . Значения  также приведены в таблицах. Согласно таблицам первый и второй способы выбора значений  обеспечивают точность определения массы в несколько процентов. Следует отметить, что использование при вычислении  полинома первой степени дает практически те же значения , что и пятой степени, но с меньшей погрешностью . Так в случае полинома первой степени для способов 1 и 2 имеем = и = соответственно; в случае полинома пятой степени для тех же способов = и = (см. табл. 2). Использование полиномов первой и пятой степени при оценке  дает практически одинаковый результат – =для способа 1 и = для способа 2.

Третий и четвертый способы выбора  оказались заметно точнее. Погрешность  способа 4 при значениях  Н и  Н оказалась всего 0.16 % и 0.13 %, погрешность  не превысила 1 %.

Результаты оценок массы при использовании значения , полученного в ходе анализа телеметрической информации о давлениях в магистралях ДУ, представлены в табл. 4. Описанные выше способы 1, 2, 3 и 4 оценивания массы применялись с единственным изменением: в формулах (2.1), (2.2) учитывались лишь те измерения, для которых  с. В частности, в способах 3 и 4 учитывались данные, относящиеся к пяти из десяти упоминавшихся выше интервалов постоянства функции . В данном случае оценки массы станции получились несколько точнее. В частности способ 3 дал погрешность 0.5%. К сожалению, несовершенство такого способа оценки  заключается в принципиальной невозможности уменьшить . Минимальное значение , обусловленное конструкцией двигательной установки, составляет не менее 2 %. В связи с этим, оценка  в данном случае значительно хуже соответствующих оценок, приведенных в табл. 2, 3. Однако оценки  в табл. 4, наоборот, значительно лучше аналогичных оценок в табл. 2, 3.

Таким образом, оценки значения массы при использовании значения , полученном сглаживанием данных БИПС полиномом пятой степени (182120 кг) или сообщенном двигателистами (182376 кг), оказались достаточно близкими к проектной оценке массы станции (183311 кг).

5. Исследование низкочастотных колебаний корпуса станции, вызванных работой двигателей ДПО. Для исследования были взяты шесть отрезков данных измерений MAMS, полученных 20 и 26 августа 2004 г. Два отрезка получены во время коррекций орбиты в указанные дни, остальные отрезки – при отсутствии динамических операций и включений двигателей управления ориентацией. Выбранные отрезки данных представлены на рис. 2 – 4. Рис. 2 и 3 иллюстрируют отрезки 1 и 2, на которых выполнялась коррекция орбита. Они устроены следующим образом. Рисунки 2а и 3а аналогичны рис. 1б. На них приведены графики данных измерений (пунктирные линии) и графики соответствующих аппроксимирующих выражений  (сплошные линии). Для первого отрезка , , для второго отрезка , . На 2б и 3б приведены графики разностей  – ломаные линии с вершинами в точках . Поскольку при выполнении коррекции тяга двигателей устанавливается не сразу, начальные участки рядов данных  были исключены из анализа. В случае отрезка 1 были отброшены данные, отвечающие первым 30 с, в случае отрезка 2 – данные первых 20 с. В обоих случаях для анализа были выбраны 2048 последовательных значений разностей . Графики подвергнутых анализу значений  приведены на рис. 2в и 3в.

На отрезках 3 – 6 включений двигателей не было, поэтому данные этих отрезков, содержащих каждый по 2048 значений (рис. 4), анализировались непосредственно. Отрезки 3, 4 расположены вблизи отрезка 1, отрезки 5, 6 – вблизи отрезка 2.

Спектральный анализ данных на выбранных интервалах выполнялся по следующей схеме. Для определенности будем говорить о данных отрезков 3 – 6. Пусть ,   – данные измерений микроускорения. Попытаемся аппроксимировать функцию  при  функцией

                                  (5.1)

где , ,  и  – параметры. Значения этих параметров будем искать методом наименьших квадратов. Составим выражение

.                                             (5.2)

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

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

,                   (5.3)

где ,   – постоянные параметры. Значения параметров находятся методом наименьших квадратов – из условия минимума функции, заданной соотношениями (5.2) и (5.3).

Чтобы проверить найденное решение другим методом, наряду с функцией  рассматривалась функция

,   ,

называемая периодограммой Шустера [5]. Пусть исследуемая функция  имеет вид (5.3) при . Тогда , периодограмма имеет локальные максимумы в точках , причем  . Точность выписанных соотношений увеличивается с ростом . Таким образом, знание максимумов периодограммы позволяет получить оценки частот и амплитуд гармонических составляющих функции .

Ниже для удобства вместо графиков функций  и  приводятся графики функций

,     .

Минимумы функции  дают оценки среднеквадратической ошибки аппроксимации функции  выражением (5.1), максимумы функции  – оценки амплитуды .

На рис. 5 приведены графики функций  и  для отрезков 1, 3 и 4, на рис. 6 – аналогичные графики для отрезков 2, 5 и 6. Точки значимых минимумов функции  и точки значимых максимумов функции  отличаются не более, чем на несколько единиц  Гц. Судя по рисункам, каждая из функций  или  содержит значительное число гармоник. Ниже в соответствующих выражениях (5.3) будем учитывать все наиболее значимые гармоники (отвечающие экстремумам функций ,  с резко выделяющимися экстремальными значениями) и несколько менее значимых. Амплитудой гармоники с частотой  в (21) будем называть величину . В качестве единицы измерения функций  и  примем  м/с, частоты будем указывать в герцах.

Частоты и амплитуды некоторых наиболее значимых гармоник, обнаруженных в функции  для отрезков 1 – 6 приведены в таблице 5. Таблица устроена следующим образом. Ее верхняя строка содержит номера исследованных отрезков данных. Каждому отрезку отвечают два столбца – столбец частот () и столбец амплитуд (). В этих столбцах указаны частоты и амплитуды построенного для функции  данного отрезка данных аппроксимирующего выражения (5.3). Параметры гармоник, обнаруженных в разных отрезках данных, но имеющих примерно одинаковые частоты, помещены в одну и ту же горизонтальную строку. Такая организация таблицы, во-первых, удобна для анализа. Во-вторых, в рассматриваемом случае естественно предположить, что примерно одинаковые частоты, проявляющиеся в разных отрезках, – результат воздействия одной и той же частоты некоторого реального фактора. Небольшие вариации частот, превышающие погрешности, обусловленные ограниченностью длины отрезков данных, могут быть объяснены нестабильностью функционирования бортового оборудования.

Номинально акселерометр MAMS позволяет определять микроускорения в диапазоне частот от 0 до 1 Гц. Однако головки прибора имеют скорость выборки 10 измерений в секунду, поэтому неудивительно, что в спектре данных измерений присутствуют частоты до 5 Гц. Сколько-нибудь значимые возмущения наблюдаются только до 1.6 Гц, поэтому периодограммы представлены на рисунках в диапазоне от 0 до 1.6 Гц. Необходимо отметить, что поиск гармонических составляющих в данных измерений является коварной задачей. Гармоники с малыми амплитудами могут порождаться случайными ошибками в данных [4]. Для гармоник с большими амплитудами таких сомнений не возникает. В рассматриваемой задаче можно быть достаточно уверенными в существовании гармоник с амплитудами более  м/с.

Некоторые группы экстремумов функций  и  с близкими частотами (см. рис. 5, 6) отвечают колебаниям с так называемым непрерывным спектром [4]. Форма представления таких колебаний в виде (5.3) не адекватна. Однако, если указанные группы экстремумов достаточно локализованы, отвечающую им составляющую функции  можно приближенно представить одной гармоникой. Именно такой подход был использован при составлении табл. 5.

Аналогичным образом проводился спектральный анализ данных измерений MAMS компонент кажущегося ускорения станции вдоль осей  и  строительной системы координат СМ. Эти компоненты обозначим соответственно  и . Введенные выше функции ,  для данных измерений величин  и  обозначим ,  и , . Графики функций , ,  и , ,  для отрезков 1, 2, …, 6 изображены на рис. 7 – 10 и 11 – 14. Частоты и амплитуды наиболее значимых найденных гармоник приведены в табл. 6, 7. Данные измерений  и  на отрезках 1, 2 (см. верхние графики на рис. 7, 11) содержат значительные выбросы, обусловленные включением ДПО. Перед проведением спектрального анализа эти выбросы были заменены средними значениями соответствующих данных на отрезке (см. нижние графики на рис. 7, 11).

Анализ отрезков 1, 2, 3, 4 и 6 позволил установить следующее. При коррекции орбиты станции с помощью ДПО грузового корабля станция испытывает значительные колебания с частотами 0.022÷0.056 Гц. При выключенных ДПО колебания с этими частотами практически отсутствуют.

Частоты, превышающие 0.6 Гц, в спектрах отрезков данных 1, 2, 3, 4 и 6 практически одинаковы. Амплитуды соответствующих гармоник на отрезках 1, 2 (во время работы ДПО) существенно выше, чем на отрезках 3, 4 и 6. Как видно из табл. 5, 6, 7 и рис. 5, 6, 9, 10, 13 и 14, амплитуды вибрационных микроускорений во время работы ДПО в среднем возрастают на порядок, а амплитуды гармоник -компоненты микроускорения с частотами 0.100÷0.113, 0.252, 0.274, 0.352 и 1.017÷1.043 Гц возрастают на два порядка. В случае -компоненты при работе ДПО наиболее значительные возмущения наблюдаются на частотах 0.101÷0.111, 0.254, 0.275, 0.360, 0.401 и 0.701÷0.719 Гц. В случае -компоненты – на частотах 0.250, 0.381÷0.425, 0.450, 0.652, 0.677, 0.805 Гц, причем в -компоненте амплитуды возрастают на два порядка.

Отрезок 5 следует рассмотреть особо. На этом отрезке происходила передача управления станцией от американского сегмента к российскому, что вызвало некоторое изменение работы системы управления ориентацией и несколько повысило уровень микроускорений. Амплитуды гармонических составляющих в данных отрезка 5 в несколько раза больше соответствующих амплитуд на отрезках 3, 4 и 6. В отличие от отрезков 1, 2, 3, 4 и 6 в спектре -компоненты микроускорения на отрезке 5 нет возмущений с частотами более 1.03 Гц, а в спектрах - и -компонент не наблюдаются возмущения с частотами более 0.808 Гц. На отрезке 5 в спектре -компоненты присутствуют значительные возмущения на частоте 0.385 Гц, в спектре -компоненты – на частотах 0.385, 0.390, 0.738 Гц, в спектре -компоненты – на частотах 0.751, 0.759 Гц.

         Данная работа выполнена при финансовой поддержке РФФИ (проект 05-01-00451).

Литература

1.     Семенов Ю.П., Григорьев Ю.И., Соловьев В.А., Благов В.Д., Курилов В.А., Бодин Н.Б., Ковтун B.C., Почукаев В.Н., Поляков B.C. Способ определения массы сборки космических объектов в процессе изменения параметров орбиты. Авторское свидетельство RU 2098326 С1.

2.     Линник Ю.В. Метод наименьших квадратов и основы математико-статистической теории обработки наблюдений. М.: ГИФМЛ, 1958.

3.     Бабкин Е.В., Беляев М.Ю., Ефимов Н.И., Сазонов В.В., Стажков В.М. Определение квазистатической компоненты микроускорения, возникающего на борту Международной космической станции. Космические исследования, 2004, т. 42, № 2, с. 162 – 171.

4.     Теребиж В.Ю. Анализ временных рядов в астрофизике. М.: Наука, 1992.