Известия РАН. Теория и системы управления, 2023, № 3, стр. 106-119

Введение. В последние годы малоразмерные беспилотные летательные аппараты (ЛА) получили значительный стимул к развитию, в том числе в связи с ростом вычислительных мощностей микроконтроллеров [1] и систем навигации [2].

Учитывая имеющиеся габарито-массовые характеристики таких аппаратов, накладываются ограничения на установку датчиков и возможность их резервирования [3]. В этой связи перспективным подходом является комплексная обработка сигналов. В общем случае она позволяет уменьшить количество необходимых на борту средств измерения, используя информацию, которая уже имеется на борту, и возросший потенциал вычислителей [4, 5]. Комплексирование также позволяет восстанавливать измерительную информацию в случае отказа датчиков, а также при расследовании летных происшествий.

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

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

1. Общий подход. Для нахождения значений углов ориентации применяется подход, использующий идеи прямого метода решения задач поиска оптимального управления [6–10]. В его рамках искомые сигналы задаются конечным числом параметров, что переводит решаемую задачу в класс задач численной оптимизации.

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

В рассматриваемой задаче предполагается, что значения проекций скорости в нормальной земной системе координат и значения перегрузок в связанной измеряются соответственно спутниковой навигационной системой и бортовыми акселерометрами, т.е. считаются заданными. Согласно известным уравнениям динамики полета [11], они связаны между собой однозначными математическими соотношениями, которые определяются значениями углов ориентации. Это предоставляет возможность получить оценки углов ориентации.

Так как в общем случае искомые сигналы могут определяться большим числом параметров, то для решения поставленной задачи целесообразно использовать популяционные или генетические алгоритмы [12–14]. В статье применяется метод роя частиц [6, 15].

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

Решение задачи основывается на известных формулах, связывающих между собой нормальную земную и связанную системы координат [16].

Пусть известны значения перегрузок в связанной системе координат. Тогда можно определить проекции ускорения на оси этой системы [11]:

где ${{a}_{{xg}}},\;{{a}_{{yg}}},\;{{a}_{{zg}}}$ – значения проекций ускорений, приведенные в нормальной земной системе координат; ${{n}_{x}},\;{{n}_{y}},\;{{n}_{z}}$ – перегрузки, измеренные относительно осей связанной системы координат; $\vartheta ,\;\gamma $ – углы тангажа и крена.

Воспользовавшись матрицей перехода из связанной системы в нормальную земную [16], можно пересчитать ускорения в этой системе:

где $\psi $– угол рыскания. Проинтегрировав их, получим значения проекций скоростей: где ${{V}_{{xg}}},\;{{V}_{{yg}}},\;{{V}_{{zg}}}$ – значения проекции вектора скорости в нормальной земной системе координат; t₀ – время начала интервала.

Эти операции можно производить и в обратном порядке: выполнить интегрирование скоростей в связанной системе координат:

где ${{\omega }_{x}},{{\omega }_{у}},{{\omega }_{z}}$ – угловые скорости относительно осей связанной системы координат, и перевести скорости в нормальную земную систему:

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

В основе подхода лежит идея, что измерения перегрузок и скоростей могут быть получены с высокой точностью и от не зависимых друг от друга систем, использующих различные физические принципы: ускорения определяются по показаниям акселерометров, а проекции скоростей в земной системе – по данным спутниковой системы навигации. В то же время формулы пересчета определяются значениями углов ориентации ЛА, прежде всего тангажом и креном. Это позволяет предположить, что, добившись хорошего соответствия между перегрузками и проекциями скорости, удастся определить значения Эйлеровых углов, дифференцируя которые можно получить оценки угловых скоростей [11].

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

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

где α, β – углы атаки и скольжения, $\Theta ,\;\Psi $ – углы наклона и поворота траектории.

Из равенств (1.6), (1.7) видно, что они зависят от значений угла рыскания. Кроме того, они связывают значения углов атаки и скольжения, которые задают направление вектора скорости относительно связанной системы координат, с углами наклона и поворота траектории $\Theta $ и $\Psi $, определяющими направление вектора скорости в земной нормальной системе координат [16].

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

Более подробное описание можно найти в [18]. Аналогичные формулы могут быть применены и для оценки углов наклона и поворота траектории, если подставить в них значения проекций скорости в земной системе координат.

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

В общем виде полученная оптимизационная задача запишется так:

где ${\mathbf{x}}$ – вектор фазовых переменных системы, определяемой уравнениями (1.1)–(1.3) или (1.4), (1.5); ${\mathbf{u}}$ – вектор искомых сигналов, аналогичный вектору управления в классических задачах оптимального управления; $U$ – множество возможных значений управляющих сигналов.

Для описанной выше конкретной постановки функционал принимает вид

где ${{t}_{0}}$, $T$ – время начала и конца участка соответственно; ${\mathbf{u}} = {{\left[ {\vartheta \;\gamma \;\psi } \right]}^{{\;{\text{T}}}}}$ – вектор углов ориентации; ${{{\mathbf{V}}}_{{\mathbf{g}}}} = {{\left[ {{{V}_{{xg}}}\;{{V}_{{yg}}}\;{{V}_{{zg}}}} \right]}^{{\;{\text{T}}}}}$ – вектор скорости в земной нормальной системе координат, ${{{\mathbf{\tilde {V}}}}_{{\mathbf{g}}}}$ – известный вектор скорости в земной нормальной системе координат.

Для формулы (1.11) U представляет собой множество возможных значений углов ориентации ϑγψ. Границы этого множества задаются априорно при помощи системы неравенств по каждому углу.

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

Будем предполагать, что при этом угол рыскания является известным. Обозначим такую формулировку задачи (1.12).

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

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

В данном исследовании были применены кубические сплайны дефекта 2, которые также называются Эрмитовыми сплайнами. Их особенность заключается в разрывности второй производной, поэтому требуется задавать в узлах сплайна значения как самой функции, так и ее первой производной. Более подробно этот вид сплайнов описывается в [19].

Данная параметризация не является единственно возможной. Можно, например, предложить использовать интерполяционные кубические сплайны [19], имеющие в узлах сплайна непрерывную вторую производную. Это условие задает систему уравнений, которая разрешается относительно значений производной в узлах. Она описывается трехдиагональной матрицей, и для ее решения существует хорошо разработанный способ, называемый “методом прогонки” [19].

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

Сравнение этих двух видов сплайнов для рассматриваемой задачи проводится в разд. 2.4.

1.3. Метод роя частиц. Предлагаемый метод преобразует задачу нахождения сигнала к задаче численной оптимизации. Для ее решения предлагается использовать алгоритм, называемый методом роя частиц, идея которого впервые была предложена в [20].

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

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

На настоящий момент существует множество модификаций алгоритма роя [14, 22]. Описание основных идей алгоритма приводится в [15, 23]. Реализованный в статье алгоритм подробно описан в [24].

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

Количество узлов сплайна, необходимое для описания управляющего сигнала, определялось путем решения задачи аппроксимации. При рассмотрении участков продолжительностью 7 с для задания углов тангажа и крена применялись сплайны из 8 узлов, а угла рыскания – из 4. При этом требовалось решить задачу оптимизации для 40 параметров, что персональный компьютер, оснащенный процессором Intel Core 2 Duo с частотой 2.4 ГГц и оперативной памятью 1.99 Гб, выполнял за 20–30 мин.

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

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

С другой стороны, разрывность второй производной Эрмитова сплайна также может создавать дополнительные погрешности. Как можно видеть на графике угловой скорости рыскания (второй график рис. 1), при избыточной плотности узлов это приводит к возникновению заметных искажений сигнала (интервал 4.5–5 с).

2.2. Выбор математической модели. Чтобы сравнить модели (1.1)–(1.3) и (1.4), (1.5), рассматривалась задача с известным значением угла рыскания, т.е. в постановке (1.12). Результаты сравнения приводятся на рис. 2.

Можно видеть, что интегрирование скоростей по уравнениям (1.4), (1.5) дает несколько худшие результаты по сравнению с интегрированием скоростей в нормальной системе. Это подтверждается и значениями среднеквадратичных отклонений (СКО), приведенных в табл. 1. Поэтому в дальнейшем приводилось решение задачи для модели объекта, задаваемой уравнениями (1.1)–(1.3).

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

2.3. Восстановление значений углов ориентации по информации о проекциях перегрузки и скорости. При решении задачи поиска трех углов ориентации ЛА было установлено, что в рамках рассматриваемой модели не удается однозначно их восстановить. Оценки тангажа, крена и рыскания, полученные с использованием соотношений (1.1)–(1.3), содержат существенные погрешности, особенно по углу рыскания. Это иллюстрируется рис. 3.

Как предполагалось выше, сформулированная на основе модели (1.1)–(1.3) и целевого функционала (1.11) задача получается плохо обусловленной, поэтому решение не удается восстановить однозначно. Это приводит в первую очередь к существенным погрешностям при определении сигнала рыскания (рис. 3).

Рассмотрим способы повышения качества оценок углов, привлекая для этого дополнительную информацию. Воспользуемся выражением (1.7). Для этого добавим в состав математической модели оценки углов атаки и скольжения по формулам (1.8) и (1.9), а также значения углов наклона и поворота траектории по аналогичным формулам. Это позволяет оценить разность между левой и правой частями равенства (1.7) и внести ее в функционал с некоторым весовым коэффициентом. Целевой функционал приобретает вид

где k₂ – весовой коэффициент, подбираемый эвристически,

Как следует из рис. 4, использование функционала (2.1) позволяет получить значения углов тангажа и крена достаточно близкие к эталонным, но погрешность определения угла рыскания сохраняется недопустимо высокой. Значения найденных при этом СКО приводятся в табл. 2.

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

где k₁, k₂ – весовые коэффициенты, подбираемые эвристически,

Полученные для целевого функционала (2.2) значения углов ориентации приводятся на рис. 5, соответствующие значения СКО – в табл. 2.

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

где ${{k}_{1}},\;{{k}_{2}},\;{{k}_{3}}$ – весовые коэффициенты, подбираемые эвристически, ${\mathbf{n}} = {{[{{n}_{x}}\;{{n}_{y}}\;{{n}_{z}}]}^{{\text{T}}}}$ – вектор перегрузок в связанной системе координат.

Значения весовых коэффициентов подбирались из соображений примерного равенства вклада различных слагаемых в функционал. В рассматриваемой задаче они составили k₁ = 64, ${{k}_{2}} = 64,\,\,{{k}_{3}} = 20$. Коэффициенты получаются безразмерными, поскольку ставились перед выражениями, составленными из тригонометрических функций, или безразмерными величинами.

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

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

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

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

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

2.4. Сравнение видов сплайнов. Для рассматриваемой задачи проводилось, кроме того, сравнение Эрмитовых кубических сплайнов и интерполяционных кубических сплайнов. Оба вида сплайна применялись в решении задачи при предположении об известном значении угла скольжения. Соответствующие значения СКО приводятся в табл. 4, а изображения результатов – на рис. 7.

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

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

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