Спосіб прогнозування вектора стану центра мас космічного апарата на орбіті

Реферат: Винахід належить до космічної техніки і може бути використаний для прогнозування вектора стану центра мас космічного апарата. Спосіб прогнозування вектора стану центра мас космічного апарата на орбіті, який полягає у тому, що система диференціальних рівнянь руху центра мас космічного апарата вирішують методом чисельного інтегрування Рунге-Кута. Інтегрування проводиться з кроком більшим, ніж потрібний крок формування вектора стану КА, а для його визначення з потрібним кроком застосовується екстраполяція поліноміальною функцією, яка основана на розкладанні в степеневий ряд Тейлора-Маклорена і являє собою поліноміальні залежності для радіус-вектора і вектора швидкості в кінці кроку екстраполяції так, що радіус-вектор і вектор швидкості в кінці кроку екстраполяції залежать від радіус-вектора і його похідних на початку цього кроку. Порядок поліноміальної функції відповідає максимальному порядку врахованих в ній похідних радіус-вектора. Винахід дозволяє зменшити час, що витрачається бортовою цифровою обчислювальною машиною при визначенні вектора стану центра мас космічного апарата з потрібним кроком формування, без зниження необхідної точності його розрахунку. UA 104485 C2 (12) UA 104485 C2 UA 104485 C2 5 10 15 20 25 30 35 40 45 50 55 60 Винахід належить до космічної техніки і може бути використаний для прогнозування радіусвектора і вектора швидкості центра мас космічного апарата (КА) на орбіті. Спосіб прогнозування вектора стану центра мас є важливим елементом системи управління КА. Від нього залежить точність визначення радіус-вектора і вектора швидкості, а також час, що витрачається бортовою цифровою обчислювальною машиною (БЦОМ) на визначення цих навігаційних параметрів. У свою чергу від точності, з якою розраховуються радіус-вектор і вектор швидкості центра мас КА, залежить точність побудови орієнтації і стабілізації при виконанні КА таких функцій, як спостереження за планетами і зірками, зондування земної поверхні, корекція орбіти. Тобто від вибору способу прогнозування вектора стану КА істотно залежить як технічна, так і економічна ефективність КА в цілому. Відомий ряд способів прогнозування вектора стану центра мас КА аналітичними методами. При розробці таких аналітичних методів використовують різноманітні методи наближеного інтегрування рівнянь в оскулюючих кеплерових елементах або намагаються знайти такий вид потенційної функції, що апроксимує гравітаційне поле притягаючого центра, яка допускає вирішення диференціальних рівнянь в квадратурах. Для навколоземних орбіт часто використовують способи прогнозування збуреною орбітою, коли враховується вплив другої зональної гармоніки в розкладанні потенціалу гравітаційного поля Землі в ряд по сферичних функціях [див. Дубошин Г.Н. Справочное руководство по небесной механике и астродинамике. М.: Наука, 1976, C. 563-577, Нариманов Г.С. Основы теории полета космических аппаратов. - М: Машиностроение, 1972, с. 347-349, Тихонравов М.К. Основы теории полета и элементы проектирования искусственных спутников Земли. - М.: Машиностроение, 1967, с. 109-119]. Недоліком аналітичних методів є те, що вони використовують наближені формули і не враховують широкий спектр діючих на КА збурень. А також те, що для будь-якого типу орбіти практично неможливо розробити узагальнений аналітичний метод рішення. Відомий ряд способів прогнозування вектора стану центра мас КА наближеним інтегруванням системи диференціальних рівнянь збуреного руху КА в оскулюючих кеплерових елементах орбіти. Для опису збуреного руху найчастіше використовують рівняння Лагранжа або рівняння Ньютона (див. Дубошин Г.Н. Справочное руководство по небесной механике и астродинамике. - М.: Наука, 1971, с. 266-271). Але ці системи рівнянь нелінійні як по аргументу часу, так і по шуканих функціях - оскулюючих елементах, і тому не інтегруються в кінцевому вигляді. Оскільки оскулюючі елементи є функціями, які поволі змінюються за часом, то рівняння руху в оскулюючих елементах зазвичай інтегрують наближеними методами, що є істотним недоліком прогнозування в оскулюючих кеплерових елементах, оскільки необхідна точність обчислення вектора стану КА може забезпечуватися на коротких інтервалах часу. Відомий ряд способів прогнозування вектора стану центра мас КА чисельним інтегруванням системи диференціальних рівнянь руху центра мас КА в прямокутних координатах. Серед них найчастіше використовують метод Адамса зі змінним кроком інтегрування (див. Дубошин Г.Н. Справочное руководство по небесной механике и астродинамике. - М.: Наука, 1971, с. 478-481, Нариманов Г.С. Основы теории полета космических аппаратов. - М.: Машиностроение, 1972, с. 333-344). Недоліком способів прогнозування чисельним інтегруванням в прямокутних координатах є залежність точності і часу розрахунку вектора стану КА від повноти обліку збурюючих чинників в правих частинах диференціальних рівнянь і від величини кроку інтегрування, а також те, що для забезпечення високої точності розрахунку вектора стану необхідно використовувати відносно невеликі кроки інтегрування. Відомий найбільш близький по технічній суті спосіб прогнозування вектора стану центра мас КА, який полягає в тому, що систему диференціальних рівнянь руху центра мас КА вирішують методом чисельного інтегрування Рунге-Кута з постійним кроком в прямокутних координатах. Особливістю такого методу є те, що він дозволяє врахувати широкий спектр діючих збурень, а також те, що витрати часу бортового обчислювача на виконання інтегрування пропорційні кількості реалізованих кроків інтегрування. Недоліком цього способу є те, що використання його з малим кроком формування вектора стану центра мас КА призводить до високих витрат часу БЦОМ, необхідного для виконання розрахунків. Задачею запропонованого винаходу є істотне зменшення часу, що витрачається БЦОМ на обчислення вектора стану КА способом прогнозування чисельним інтегруванням методом Рунге-Кута, що дозволяє заощадити ресурси БЦОМ для вирішення інших завдань. При цьому забезпечується необхідна точність розрахунку радіус-вектора і вектора швидкості КА. Поставлена задача вирішується тим, що чисельне інтегрування методом Рунге-Кута проводиться з кроком більшим, ніж потрібний крок формування вектора стану КА, а для його визначення з потрібним кроком застосовується екстраполяція поліноміальною функцією. 1 UA 104485 C2 5 10 15 20 25 30 35 40 Поліноміальна функція основана на розкладанні в степеневий ряд Тейлора-Маклорена і являє собою поліноміальні залежності для радіус-вектора і вектора швидкості в кінці кроку екстраполяції. Радіус-вектор і вектор швидкості в кінці кроку екстраполяції залежать від радіусвектора і його похідних на початку цього кроку. Порядок поліноміальної функції відповідає максимальному порядку врахованих в ній похідних радіус-вектора. Надалі, крок, з яким використовується поліноміальна функція, називатимемо кроком екстраполяції, а крок, з яким проводиться чисельне інтегрування, - кроком інтегрування. Запропонований спосіб прогнозування можна назвати комбінованим способом прогнозування чисельним інтегруванням методом Рунге-Кута і екстраполяцією поліноміальною функцією. Як відомо, крок інтегрування може бути збільшений тільки так, щоб точність формування радіус-вектора і вектора швидкості КА впродовж необхідного інтервалу часу задовольняла встановленим вимогам. Крок екстраполяції і порядок поліноміальної функції необхідно вибирати такими, щоб забезпечувати визначення вектора стану з необхідною точністю і кроком формування. Повна похибка запропонованого способу в основному залежить від порядку методу чисельного інтегрування, кроку інтегрування і повноти обліку збурень в правих частинах диференціальних рівнянь руху центра мас КА, а на інтервалах усередині кроку інтегрування від величини похибки чисельного інтегрування, порядку поліноміальної функції, кроку екстраполяції, форми представлення похідних радіус-вектора і вектора швидкості КА в поліноміальній функції і складу врахованих збурень в її моделі прискорення і моделях похідних прискорення. Важливою перевагою такого підходу є відсутність на інтервалах усередині кроку інтегрування багатократного повторення ітераційних операцій, а також необхідність наявності тільки одного відомого вектора стану КА на початку цих інтервалів. Об'єм розрахунків радіусвектора і вектора швидкості при використанні поліноміальної функції залежить від тих же чинників, що і точність визначення вектора стану усередині кроку інтегрування. Запропонований спосіб забезпечує можливість зменшення часу, що витрачається БЦОМ на обчислення вектора стану КА з потрібним кроком формування, без зниження необхідної точності його визначення. В рамках поставленої задачі важливими характеристиками запропонованого комбінованого способу прогнозування є відносні витрати процесорного часу, які розраховуються як відношення витрат процесорного часу комбінованого способу до витрат процесорного часу якого-небудь способу прогнозування вектора стану центра мас КА. За допомогою цієї характеристики можна зробити висновок про доцільність заміни цього способу прогнозування на комбінований з погляду зменшення часу, що витрачається БЦОМ на обчислення вектора стану КА. Якщо відоме відношення витрат процесорного часу на один крок екстраполяції поліноміальною функцією до витрат на один крок прогнозування яким-небудь іншим способом ΔРТе, і відоме таке ж відношення для інтегрування методом Рунге-Кута APT, то відносні витрати часу комбінованого способу (при рівності кроків екстраполяції і прогнозування) розраховуються за формулою:  h  he PT i  PT e 1  e  , hi>he,  hi hi    де he - крок екстраполяції комбінованим способом; hi - крок інтегрування комбінованим способом. Іншими важливими характеристиками запропонованого комбінованого способу прогнозування є сферичні похибки відхилень по радіус-вектору ΔR і вектору швидкості ΔV, які показують, наскільки знижується точність розрахунків цих навігаційних параметрів КА на ділянках використання поліноміальної функції в порівнянні з еталонними значеннями. Як початкові умови для формування еталонних значень використовуються радіус-вектор і вектор швидкості КА, що проінтегровані на початок інтервалів усередині кроку інтегрування. Це дозволяє виключити з аналізу вплив похибок задання початкових умов на точність екстраполяції поліноміальною функцією. Метод, вибраний для формування еталонів, необхідно використовувати з кроком меншим або рівним кроку екстраполяції. Сферичні похибки відхилень розраховуються за формулами: PT  45 50 R  55 xк  х ет 2  ук  у ет 2  zк  z ет 2 V  v хк  v хет 2  v у к  v у ет 2  v zк  v z ет 2 ; , де хет, ует, zет - компоненти еталонного радіус-вектора; хк, ук, zк - компоненти радіус-вектора, розраховані комбінованим способом прогнозування; v хет, vyeт, vzeт - компоненти еталонного 2 UA 104485 C2 5 10 15 радіус-вектора; vхк, vук, vzк - компоненти вектора швидкості, розраховані комбінованим способом прогнозування. Суть винаходу пояснюється кресленнями, де зображені: на фіг. 1 - характеристики комбінованого способу прогнозування чисельним інтегруванням методом Рунге-Кута четвертого порядку і екстраполяцією поліноміальною функцією другого порядку, на фіг. 2 - характеристики комбінованого способу прогнозування тим же чисельним методом інтегрування і екстраполяцією поліноміальною функцією третього порядку, на фіг. 3 - характеристики комбінованого способу прогнозування тим же чисельним методом інтегрування і екстраполяцією поліноміальною функцією четвертого порядку. Для демонстрації працездатності та ефективності запропонованого способу прогнозування вектора стану центра мас розглянемо три його різновиди з різними кроками екстраполяції і інтегрування. Розрахунок вектора стану проводиться в Грінвічській системі координат. У всіх цих прикладах використовується метод Рунге-Кута четвертого порядку. У системі диференціальних рівнянь руху КА враховується вплив в розкладанні потенціалу гравітаційного поля Землі U в ряд по сферичних функціях гармонік до шістнадцятого порядку включно (див. Астрономический ежегодник на 2000 год. Институт прикладной астрономии РАН, 1999):  16  a  n n  k 1    r   Cnk cos k  Snk sin k Pn sin  ,  r  n2   k 0   де μ - геоцентрична гравітаційна постійна Землі; α - велика піввісь загального земного еліпсоїда; r - модуль радіус-вектора КА; Cnk, Snk - коефіцієнти гармонійного розкладання; λ, φ Ur, ,   20 25 30 35 40 k довгота і широта підсупутникової точки; Pn - приєднані функції Лежандра. Розглянуті приклади відрізняються порядком поліноміальних функцій: другий, третій і четвертий. Для поліноміальної функції другого порядку застосовувалася наступна формульна схема:    v t n  t   v t n   gt n t;      t r t n  t   r t n   v t n   v t n  t  , 2     де, r t n  , v t n  і gt n  - радіус-вектор, вектор швидкості і вектор прискорення у момент часу tn; Δt - крок екстраполяції. Формульна схема для поліноміальної функції третього порядку:    t    v t n  t   v t n    gt n   gt n  t; 2        t  t    t    v t    gt   gt  t  t t, r n r n  n  n  n     3 2       t  - похідна вектора прискорення у момент часу tn. де g n Формульна схема для поліноміальної функції третього порядку:     t  t     v t n  t   v t n    gt n    gt n   gt n   t;  3 2           t  t    t    v t    gt    gt   t  t  t  t t. r n  n  n  n g n     r n  4 3 2        де t n  - друга похідна вектора прискорення у момент часу t n. g У моделях прискорення і моделях похідних прискорення, використаних в поліноміальних функціях, враховується вплив в розкладанні потенціалу гравітаційного поля Землі тільки другої зональної гармоніки. Числові розрахунки проводяться для орбіти, параметри якої протягом доби змінюються так: модуль радіус-вектора від 7022,966 км до 7081,509 км, модуль вектора швидкості від 7,489 км/с до 7,551 км/с, ексцентриситет від 0,003191 до 0,005188, нахил орбіти від 98,2654° до 98,2806°, довгота висхідного вузла від 120,0302° до 121,0817°. Проаналізуємо відносні витрати процесорного часу розглянутих різновидів комбінованого способу. Як еталонний спосіб прогнозування виберемо такий же метод чисельного інтегрування, як і в даних прикладах. При цьому ΔРТi=1, а ΔРТe=0,129 для поліноміальної функції другого порядку, ΔРТе=0,266 для поліноміальної функції третього порядку, ΔРТe=0,475 3 UA 104485 C2 5 10 15 20 25 30 для поліноміальної функції четвертого порядку. Вибираючи різні кроки екстраполяції і інтегрування, отримуємо різні значення відносних витрат часу (фіг. 1-3). У всіх випадках ΔРТ