Определение статистических характеристик движения самолета при автоматическом заходе на посадку Текст научной статьи по специальности «Физика»

УЧЕНЫЕ ЗАПИСКИ НАГИ Том XVI 1985

№ 6

УДК 629.735.33.051.83—52

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

В. П. Кузьмин, Г. В. Парышева-

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

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

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

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

л

нии [1]. Измеренное угловое отклонение е определяется выражением

(/>)[*+/,(*)], (1)

где №пр(р) — передаточная функция приемника, р — переменная Лапласа, 5 — случайный коэффициент, определяющий разброс крутизны

глиссадной или курсовой зон [1], /е (0—случайная функция времени, моделирующая радиотехнические помехи, е— истинное угловое отклонение, связанное с линейным отклонением высоты (в случае продольного движения)

(2)

где Ь — расстояние от антенны приемника на самолете до маяка.

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

/. = 0§/0(£).£(*), (3)

где — заданная функция дальности до маяка, а, — случайный

параметр, /¡(¿)— стационарный гауссовский случайный процесс.

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

их(Н) = ихФ(Н); |

иж{Н) = и,Ф(Н), | ^

где Ф(Я) =0,22 1пЯ+0,5; Я— высота центра масс самолета в метрах, Цх, — независимые случайные параметры, определяющие значения компонент ветра на высоте Я= 10 м, измеряемые в м/с. Конкретные законы распределения величин их и иг приведены ниже {см. (28)].

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

Для модели атмосферной турбулентности, рекомендованной ИКАО [2] и принятой при расчетах в данной работе, интенсивность компонент турбулентности связана с величиной модуля градиентного ветра

0,18«, 1

аГу = 0,09и, } (5)

где о\ух, —средние квадратические отклонения продольной,

боковой И вертикальной компонент соответственно, М = ]/«2 + —

случайная для каждой реализации величина модуля ветра на высоте Я = 10 м.

Таким образом, все три компоненты турбулентности могут быть представлены в виде, аналогичном (3), с той лишь разницей, что интенсивности компонент турбулентности не являются независимыми, а определяются соотношениями (5).

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

7.(0, ^*(0, %{?), ^(0-

Рассматривается задача определения статистических характеристик для модели эквивалентной системы [3]. В этом случае система уравнений включает в себя дополнительные уравнения для формирую-

щих фильтров, возбуждаемых белыми шумами, для определения случайных функций /Е(0. Wx{t), Wy(t) и Wz(t). В дальнейшем рассматривается продольное движение самолета, так как структура уравнений бокового движения аналогична рассматриваемой.

Эквивалентная линейная система в случае продольного движения самолета с учетом соотношений (1) — (5) будет иметь вид

=A{S, t)x + V'ul + á\ В (0-6, + С (t) Ъ + D (t) п., (6)

dt

где x — вектор фазовых координат эквивалентной системы, и |2 — векторы независимых гауссовских белых шумов единичной интенсивности, Л(5, t), B(t), C(t), D(t) — матрицы соответствующих размерностей, полученные при линеаризации уравнений управляемого движения самолета.

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

+ (7)

dt

и может быть записано уравнение для корреляционной матрицы К вектора x(t) [3]:

= A (i) • К(Í) + К(t) Лт (Í) + В (t) L (t) Вт (*), (8)

dt

где L (í) — матрица интенсивностей компонентов белого шума §.

Матрица к. имеет размерность пХп, где п — размерность вектора х(1), но в силу симметрии матрицы К решение уравнения (8) эквивалентно решению системы П ,П^ ^ обыкновенных дифференциальных уравнений первого порядка.

Решение систем уравнений (6), (8) подразумевает задание начальных условий х(0) и /С(0). При решении задач на ограниченном интервале времени матрица /С(0) не известна, так как она является результатом действия возмущений в предшествующие рассматриваемым моменты времени. Неизвестным также является вектор х(О) при наличии постоянно действующего возмущения их.

Если рассматриваемый интервал времени достаточно велик и значительно превышает время затухания переходных процессов в системах (6), (8), то влияние начальных условий для таких интервалов времени будет практически несущественным. Возможность рассмотрения изолированного участка движения по глиссаде определяется именно большой протяженностью этого участка от высоты «захвата» глиссады [1] (Я=400 м) до высот # = 60-П5 м, для которых определяются статистические характеристики точности. В этом случае в начальных условиях для корреляционной матрицы эквивалентной системы следует задавать только элементы, соответствующие возмущениям; все прочие элементы можно задавать равными нулю. Для вектора х можно также принять х(0) =0.

2. Определение корреляционной матрицы. При определении статистических характеристик вектора фазовых координат x(t), описываемого уравнением (6), рассмотрим случай постоянного значения крутиз-

ны 5. Определим последовательно изолированное влияние трех слагаемых, представляющих возмущения в уравнении (6) и соответствующих воздействию турбулентности Vи2х-\-u,2z В (t)-%,, радиотехнических помех оeC(tf)-g2 и горизонтального градиентного ветра D (t) их. Корреляционные матрицы, соответствующие воздействию турбулентности или радиотехнических помех, определяются при постоянных значениях иио£ соответственно. Так, если /С10 (5) — корреляционная матрица, определенная при значении и = и0, то при произвольном значении и корреляционная матрица определяется соотношением

KAS) = K10 (5)-4-. (9)

“о

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

K2(S) = K20(S)4~- (10)

°£ 0

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

чески удобнее определять воздействие параметра их ферез функции влияния, так как для определения векторной функции h(t) необходимо интегрировать систему уравнений (6), имеющую значительно меньшую размерность, чем система уравнений (8). Пусть h(S, t) — векторная функция влияния постоянного параметра иж0 на интересующие нас фазовые координаты. Тогда значения фазовых координат определяются соотношением

x(t) = h(S, €)-$*-. (И)

их о

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

Af[*(01 = jA(S, t)^fs(S)dS,

где }s(S)—плотность распределения вероятностей крутизны, М[ ] — символ математического ожидания.

В дальнейшем рассматриваются статистические характеристики вектора x(t) в фиксированный момент времени и переменная t опускается.

Соотношения (9), (10) и (11) позволяют определить случайные значения вектора х при известных законах распределения случайных параметров S, их, их и ое:

_

x = T(S) ъ + Т2 (5) -р- h (S) , (12)

и0 °Е 0 ихо

где Г17'1 = /Сю, T2Tz = K2Q, i'll и ^2 — случайные векторы, компоненты которых — нормально распределенные независимые числа с единичной дисперсией и нулевым математическим ожиданием.

Соотношение (12) предполагает, что для каждого случайного значения 5 определяются матрицы /С10 и К2п путем решения системы (8) и определяется вектор Ь(5). Соотношение (12) позволяет определить корреляционную матрицу и плотности распределения вероятностей компонент вектора х.

При фиксированном значении 5 получим

= + (13)

“о

где Оа— дисперсия параметра их.

С учетом случайности крутизны 5 корреляционная матрица определяется формулой

К = $К (5)4(5)^. (14)

Из соотношения (13) следует, что при фиксированном 5 корреляционная матрица К(5) может быть определена в процессе одного интегрирования уравнения (8) с учетом всех возмущений при постоянных значениях случайных параметров п и аЕ, определяемых соотношениями

«о = м К1 = // /«!+“* /их (Ид) Лг (Иг) йих ёиг;

°10 = М[о1] = 1^{*г)йог, (15)

где /„х, faz м fa — плотности распределения вероятностей параметров их, иг и интенсивности помех ое соответственно.

Соотношение (12) позволяет также определить распределения вероятностей компонент вектора х. Для определения распределения вероятностей следует учесть, что при фиксированных значениях параметров 5, их, иг и условное распределение фазовых координат движения самолета определяется только гауссовскими случайными возмущениями и, следовательно, является нормальным. Для фиксированного значения 5 плотность вероятности произвольной координаты определяется СООТНОШеНИеМ

/(лг/5) = 111е 2а’ Лх(их) (°г)йихйиг(1з£, (16)

/I2 1 **2 2

_______^(5) иг | гЯ) ае

где о — дюн--------—2--[-Дго//—2—•

и0 °е0

Безусловная плотность распределения вероятностей определяется формулой

/(х1) = $/(х1/8)/';(5)с18. (17)

Формально формулы (16) и (17) для определения плотности распределения вероятностей можно объединить в одну. Практически удобнее вычисления проводить последовательно. Учитывая трудоемкость определения матриц /Сю (5) и /С20 (5), при численных расчетах непрерывное распределение величины 5 заменяется дискретным,

определяемым методом эквивалентных возмущений [3]. В этом случае определение интегралов (14) и (17) сводится к определению конечных сумм.

3. Определение законов распределения параметров траектории самолета. Рассмотрим вид законов распределения фазовых координат в некоторых частных случаях. Рассмотрим воздействие только атмосферных возмущений и пусть компоненты градиентного ветра их и их распределены нормально с нулевыми математическими ожиданиями и одинаковыми дисперсиями. В этом случае величина м—м2 распределена по закону Релея:

(18)

гДе 0 = % = V

Предположим также, что для некоторых фазовых координат изолированное влияние градиентного ветра несущественно, т. е. /гг = 0. Тогда при фиксированном значении 5 плотность распределения вероятностей фазовой координаты Х{ определяется соотношением

оо

/(■*//5)= I ---------/ 1 7x1/2 ехр//---------О9)

где Ц, = Км и ■

Дисперсия координаты Х{ с учетом формул (13) и (18) определяется соотношением

(20)

“о

Учитывая соотношение (20) и значение соответствующего табличного интеграла [5], определим окончательный вид плотности распределения вероятностей (19):

-V* Щ

/^(^/5)------4-----* ^ , (21)

У" 2

'*1

где аХ1 = V ВХ1.

Рассмотрим теперь, как изменится плотность распределения фазовых координат, если учесть еще и воздействие радиотехнических помех фиксированной интенсивности. В этом случае необходимо определить закон распределения суммы двух независимых случайных величин Хг = *1+Х2, если плотности распределения вероятностей слагаемых заданы формулами

Интегральный закон распределения вероятностей величины Хг определяется интегралом

Рассмотрим случай положительных значений у.

Введем замену переменных

— х2 — х, — у

х2 — ~, x,= — , у = —.

i °t 1 «1 ^ =i

После несложных преобразований соотношения (22) можно привести к виду

-2

+ оо х2 оо

Р(*( >})= j f + (23V

1 У~*2

где а = а21а1, а выражение для ДУ имеет вид

-2

оо _ х2 +00

AY = §^r_e 2°* J r~{e-^\^-e~^)dxtdx2. (24)

У у—¿2

Если величина у велика, то величиной AY в (23) можно пренебречь и получить приближенное соотношение

+ °° - / ~2 \ Р(х,>у)~±е-**> j ^_exp(--g--H^]d*2 =

— ОО

= -±-е-^Уе”. (25)

Нетрудно заметить, что остаточный член А У отрицательный, т. е.

формула (25) дает значение вероятности с избытком. После некоторых

преобразований выражения (24) можно получить неравенство:

00 /2 \

I I < -Г'е“ УТ} J уЩ7'0Хр (- W + V2 I dx„ (26>

У

Правая часть выражения (26) отличается от соотношения (25) только пределами интегрирования. В этом случае, если у > хч, гделгг — точка максимума подынтегральной функции в (25) и (26), то величина (АУ| по крайней мере вдвое меньше основного члена (26), так как подынтегральная функция является симметричной относительно значения х2 == х*.

Таким образом, условие для практического использования формулы (24) имеет вид:

У > (27}

На рис. 1 приведено сравнение отношения Р = Р(Х{>у)/Ро, где

Я0 = -^-ехр^—^-1/2^, определенного численным интегрированием

по формуле (23) (пунктирные линии) с приближенной формулой (25). Приведенные результаты показывают, что соотношение (27) достаточно точно определяет пределы применимости формулы (25).

4. Результаты расчетов. Рассмотрим примеры расчетов по определению точности захода на посадку пассажирского самолета для реальных моделей возмущений и алгоритмов управления. Используемые модели возмущений в основном соответствуют рекомендованным ИКАО

[2]. Случайные параметры их и uz имеют нормальный усеченный закон распределения со следующими характеристиками:

\=-2,7 м/с, Маг = 0, |

= «и, = 3,75 м/с, [ (28)

|иг|<7,7м/с; — 12,8 м/с <«_,.< 5,1 м/с. J

Под усеченным по уровням [а, Ь] распределением / (х) понимается распределение f(x), определяемое формулой

/(*)

/ (■*) =

j'/w

йх

О при х>Ь или х<С.а.

Спектральные плотности для продольной и вертикальной компонент турбулентности соответствуют модели Драйдена (6]. Радиотехнические помехи представляются двумя составляющими с фиксированными (неслучайными) интенсивностями и спектральными плотностями:

с V» (»3 + <о* + ч*)

Оо>/ -

я [<о* 4- 2о>2 (у‘ + <о|) +VJ -I- 0)/]2 62

где ¿ = 1, 2, VI = <о1 = 0,2 1/с, у2 = 1,5 1/с, ае1 = аб2г=0,03°, ш2 = = 4,7 1/с.

Функция ¡о(Ь) в (3) задается в виде

Г 1 при Ь < 1100 м,

/о (£)=■{ 1 + 8-10-5 (/. — 1100) при 1100 < Л <7300 м,

I 1,5 при £ > 7300 м.

Крутизна 5 имеет нормальный, усеченный по уровням ±3ав

[см. (29)] закон распределения с характеристиками Л1[5] = 1, сг8 = 0,16.

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

Д^ = — Ж [5] == 0; Д52,з = ±0,16; Д54,б = + 0,3;

Я, = 0,38; А = А, = 0,24; Р4 = Я5 = 0,07.

На рис. 2 приведено сравнение плотностей распределения вероят-

(закон Релея), и*=5,3м/с

I-модель ИКАО

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

Управление траекторией самолета осуществляется через стабилизацию заданного угла тангажа [4], значение которого определяется соотношением

&зад= ЧГ»(Р)г,

где величина е связана с величиной е, измеряемой приемником (1),

Л

соотношением в —К (Н) в; \Уи (Р) — передаточная функция вычислителя заданного угла тангажа.

Рассматриваются два варианта алгоритма управления, соответствующие ступенчатому и непрерывному изменению коэффициента К (//) по высоте:

[ 15 при Н > 250 м,

I вариант: К (Я) - | 6>5 при н<ЯОщ

\ 15 при //>250 м;

II вариант: К (//) = { 0 06 при я <250 м.

На рис. 3 приведены отношения средних квадратических отклонений вертикальной скорости ау , угла тангажа ап угловой скорости

У

тангажа и высоты а обусловленных изолированным воздействием

г

турбулентности и радиотехнических помех, к суммарным средним квадратическим отклонениям этих же параметров (13), (14) приведенным в табл.1

Таблица 1

м/с , М/С У я®, град ашг. град/с ая, м

I вариант 4,7 1,0 1,26 0,9 1,53

II вариант 4,7 0,52 0,36 0,3 1,49

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

В табл. 2 приведены значения коэффициентов влияния градиентного ветра на параметры движения и обусловленные этим возмущения, средние квадратические отклонения для первого варианта алгоритма и значения 5=1.

Таблица 2

V, м/с Уу, м/с град “г. град/с Н, м

Л/ 1,23 0,062 -0,009 -0.004 0,076

Ч (их) 4,65 0,23 0,034 0,015 0,28

Как следует из сравнения результатов, приведенных в табл. 1 и 2, градиентный ветер их практически полностью определяет разброс путевой скорости V и, следовательно, распределение величины будет практически таким же, как и распределение величины их. Влияние градиентного ветра на другие параметры движения практически несущественно, что подтверждает предположение, принятое при выводе приближенного соотношения (21).

На рис. 4 приведены интегральные законы распределения различных параметров движения для двух вариантов алгоритма управления, определенные с учетом всех возмущений. Характерной особенностью приведенных законов распределения параметров движения (кроме путевой скорости V) является их экспоненциальность при больших значениях отношений х = -Х‘~М ^ (см. (21)]. В этом случае, если

°х I

влияние изменения крутизны на характеристики точности невелико (II вариант) и если значения <т4 и ог определять осреднением по крутизне, приведенные результаты очень близки к получаемым по формуле (25). ,

Отметим некоторые особенности численного рещения уравнения (8) для корреляционной матрицы.

Если матрица А стационарна, то решение однородной системы

(7) имеет вид х = где Яг — собственные числа матрицы А,

I

Фг — коэффициенты. Учитывая определение корреляционной матрицы /С = М [л:-л:т], легко получить, что собственные числа матрицы правых частей системы (8) представляют собой различные суммы (Хг + ^з) собственных чисел матрицы А и, следовательно, максимальное по модулю собственное число будет вдвое больше, чем для матрицы А. Из соотношений для собственных чисел следует, что шаг интегрирования

системы уравнений (8) должен быть примерно вдвое меньшим, чем при йнтегрировании системы (7).

При исследовании реальных систем управления матрица А обычно сильно разреженная, т. е. количество ненулевых элементов матрицы А порядка п значительно меньше, чем п2. Так, в рассмотренной в данной работе задаче оценки точности продольного движения самолета при автоматическом заходе на посадку, система уравнений имеет порядок 29, однако из 841 элемента матрицы А отличны от нуля только 97. Для того чтобы исключить умножение на нулевые элементы при вычислении правых частей (8), используется специальная программирующая программа. Такая программа по маске матрицы А записывает в памяти ЭВМ процедуру умножения матриц в явном виде и автоматически опускает нулевые члены. Маска матрицы А задается двоичными константами сц, соответствующими элементами матрицы

| 1, если АуфО, сч~~\ 0, если Ац = 0.

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

ЛИТЕРАТУРА

1. Белогородский С. Л. Автоматизация управления посадкой самолета. — М.: Машиностроение, 1972.

2• Air worthiness Technical Manuel, shapter 7, section 8, part 3.

3. P о с и н М. Ф., Булыгин В. С. Статистическая динамика и теория эффективности систем управления.—М.: Машиностроение, 1981.

4. Гуськов Ю. П., 3 а г а й н о в Г. И. Управление полетом самолетов.— М.: Машиностроение, 1980.

5. Прудников А. П., Б рычков Ю. А., Маричев О. И. Интегралы и ряды. — М.: Наука, 1981.

6. Доброленский Ю. П. Динамика полета в неспокойной атмосфере. — М.: Машиностроение, 1969.

Рукопись поступила 5/VII 1984 г.