Вероятностный анализ попадания беспилотного летательного аппарата в запретную зону

УДК 629.7

А. А. ЛОБАТЫЙ, Ю. Ф. ЯЦЫНА, В. Ю. СТЕПАНОВ, А. Ю. БУМАЙ

ВЕРОЯТНОСТНЫЙ АНАЛИЗ ПОПАДАНИЯ БЕСПИЛОТНОГО ЛЕТАТЕЛЬНОГО АППАРАТА В ЗАПРЕТНУЮ ЗОНУ

Белорусский национальный технический университет

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

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

Введение

Беспилотные летательные аппараты (БЛА) приобретают всё большее распространение в различных областях человеческой деятельности, так как позволяют решать широкий спектр задач при меньших затратах по сравнению с пилотируемыми летательными аппаратами [1, 2, 3]. При этом доступность приобретения и применения БЛА для широких слоев населения создает ряд новых проблем, которые ранее не существовали. Бесконтрольное применение БЛА может представлять угрозу различным объектам, а иногда и здоровью и жизни людей. Для предотвращения этого предусмотрены различные ограничения на использование БЛА. В первую очередь определяются области воздушного пространства, в пределах которых запрещено применение посторонних БЛА. Так как проникновение БЛА в запретную зону недопустимо, то при планировании применения БЛА необходимо исключить пересечение траекторией полета БЛА границы данной зоны.

В процессе полета БЛА на его траекторию действует ряд случайных факторов: ошибки системы управления БЛА (оператора), случайные возмущения атмосферы, включающие

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

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

Для проведения исследований эволюции переменных (фазовых координат), характеризующих пространственное перемещение БЛА, необходимо иметь соответствующую математическую модель. Существует большое разнообразие форм описания математических моделей БЛА. Все они основаны на известных законах физики. Так, например, в работах [4, 5] приведены достаточно подробные математические модели БЛА, включающие как кинематические и динамические уравнения движения, так и математические выражения, описывающие работу элементов системы управления, включая уравнения для описания неопре-

деленностей и случайных воздействий. В работах [6, 7] изложены основные принципы построения компьютерных программных средств для исследования полета БЛА путём математического моделирования в пространстве состояний, то есть - моделирования всех основных элементов и подсистем БЛА.

В общем случае система уравнений, описывающая полет БЛА, как и любого летательного аппарата с учетом неопределенностей представляет собой пример функционального представления модели динамической системы как нестационарного нелинейного управляемого объекта, подверженного действию возмущений. X (0 = ф[ X ,и, X = X0,

где ф[...] - в общем случае нелинейная векторная функция, X ^) - вектор переменных пространства состояний БЛА, и^) - вектор управлений, - вектор возмущений.

Задача существенно упрощается, если уравнение движения (1) является линейным. Обоснованность линеаризации системы (1) обусловлена малой степенью отклонений параметров возмущённого движения БЛА от их номинальных значений. Обычно принимается, что если отклонения входных параметров систем 5 и ^) и вызываемые ими отклонения выходных параметров 5X ^) не превосходят 3-5% номинальных значений, то такое требование выполняется, что даёт основание представлять возмущённое движение в виде линейной комбинации и ^) = и * ^) + 5и ^), X ^) = X * ^) + 5X ^) [7]. Параметры, обозначенные звёздочкой, отвечают номинальному невозмущённому движению.

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

X (0 = А^) X ^) + В^ )и ^) + ^ Ц )£(0, X %) = X0.

Задача вероятностного анализа

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

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

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

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

Известно, что для марковского случайного процесса X ^), описываемого уравнением типа (2), плотность вероятности любой размерности полностью определяется двумя функциями - так называемой первой функцией плотности вероятности f (X, t) (для фиксированного момента времени) и условной функцией плотности вероятности f (X[X^)), которая называется плотностью вероятности перехода. Обе эти функции для многомерного (векторного) непрерывного марковского процесса X ) удовлетворяют дифференциальному уравнению Фоккера-Планка-Колмогорова (ФПК), которое для функции f (X, t) в векторно-ма-тричной форме без учёта смены структуры имеет вид [9, 10]:

д/(X,0 .

— = - ¿IV п

/ (X, = /о( Xо), где п( X, t) - вектор плотности потока вероятности вида:

п( X, t) = А( X, t) / (X, t) -- 2[У хт В( X, t) / (X, t)],

где Vх - векторный оператор дифференцирования (оператор градиента), который имеет

вид: Vx = ö// системы.

и - размерность

div n

k=1

По определению [9] вектор сноса (иногда называют вектором переноса) представляет собой скорость изменения условного математического ожидания приращения вектора фазовых координат при фиксированном значении этого вектора в момент времени t. В соответствии с определением вектор сноса вычисляется по формуле:

A

& ^0

где AX = X(t + At) - X(t) - вектор приращения фазовых координат за время &t; X - фиксированное значение вектора.

Матрица диффузии по определению есть скорость изменения дисперсии случайного марковского процесса. Поэтому её можно вычислить по формуле:

B

Заметим, что практически каждый из i-х элементов (подсистем) БЛА может быть описан стохастическим уравнением вида (2). При этом вектор состояния системы X(t) представляет собой блочную матрицу-столбец, элементами которой являются векторы состояния Xi (t) подсистем. Для процесса X (t) его полной вероятностной характеристикой является плотность вероятности f (X, t), для процесса Xi (t) - соответственно / (Xi, t).

Задача определения вероятности регистрации БЛА в запретной зоне (захвата БЛА системами обнаружения) сводится к задаче невыхода БЛА (центра масс БЛА) из запретной области из в течение времени т3 = ti+i - ti при условии нахождения его в области из в начальный момент времени t0 .

Так как в реальности время полёта и возможного нахождения в запретной зоне, как правило, существенно больше времени захвата ([tK -10] » t3), то искусственно введем в рассмотрение последовательность k = tK -10 /13 равновеликих интервалов t3 = t3i = ti+1 - ti. t0 в данном случае - момент времени, при котором БЛА попадает в область пространства из, tK - момент времени, при котором БЛА покидает область из. Для каждого момента вре-

мени ti и соответствующего интервала хзг-(/ = 0,1...£ -1) необходимо определить вероятность захвата.

Обозначим вероятность нахождения X(^ в области из в каждый текущий момент времени ti через ) = Р1(X(ti) еи3). Определим эту вероятность по формуле:

Pi(t,) = J fi(X,t)dX,

U 3

где /1(X, ti)- решение уравнения (3) при условии f (±<х>, t) = 0 .

Вероятность невыхода X ^) из области из в течение времени т3 обозначим через P2(ti + т31X(ti)еи3). В этом случае вероятность Р(^+1 - ti) захвата X^) поисковой системой с учетом вероятности срабатывания на интервале ^+1 - ^ = х3 вычисляется по формуле [8].

Р

Для определения вероятности Р2 следует проинтегрировать уравнение ФПК для плотности вероятности распределения не поглощённых реализаций /2 (X, t) с учётом их полного поглощения на границе области £/з. Воспользуемся методикой, изложенной в [11], для определения /2( X, t) из обобщенного уравнения ФПК с введением в него функций поглощения вместо границ.

fXH = -А„2( x , t) -

dt dX 2

- = /2( X, t)-Ä( X, t), P

где X, t) - нормированная функция поглощения реализаций процесса X ^), которая на границе Яи области из имеет вид:

&&,t) = 5

В выражениях (10) и (11) п2(X, t) - вектор плотности потока вероятности вида:

П2( X, t) = А( X, t) /2 (X, t) -

—B

dX

Уравнение (10) следует интегрировать при следующем начальном условии на каждом интервале [^+1 - ^ ]:

/2 (X, Ь) = -1 /1( X еи 3, ^). (13) р

Рис. 1. Сечения плотностей вероятности

Интегрируя уравнение (10) по X в бесконечной области, получим уравнение для вероятности Р2 ^) на каждом интервале [^+1 - ti ]:

Р2 ^) = —q

Здесь q(t) - интенсивность поглощения реализаций процесса X^), численно равная значению плотности потока вероятности не поглощённых реализаций, вычисленной на поглощающей границе Яи области из.

да

q

—да

Решение уравнения (14) при начальном условии Р2 (^) = 1 на каждом интервале [^ — ti ] следующее:

ti+т3

Р2(^ +Х3|X(^) е и3) = ехр

— | q(t)dt

Р3(^) = 1 — П[1 — т)Р2^, + Т IX(tг) е из)].

г=1

Таким образом, в каждый текущий момент времени ti рассматриваются две плотности вероятности: /1(X, ti) и /2(X, ti). На рис. 1 качественно показаны сечения этих плотностей при задании области из в виде заданной поглощающей границы интервала [—да, у] фазовой координаты :

При времени нахождения БЛА в запретной зоне [^ —t0]» т3 следует рассмотреть k временных интервалов, на каждом из которых события захвата считаются независимыми, а вероятность их определяется формулой (9). В этом случае вероятность захвата (вероятность невыхода БЛА из запретной области) Р3 (^) к моменту времени tr вычисляется по формуле полной вероятности, которая в данном случае имеет вид [8]:

Здесь г = 0..^ — 1; k = ^ —10 / т3.

В соответствии с формулой (17) вероятность Р3 (^) является дискретной возрастающей функцией моментов времени tr. Вероятность того, что захват произойдет к произвольному текущему моменту времени t, который не совпадает с tr, определяется интерполяцией (экстраполяцией) значений Р^г).

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

Для определения /1(X, t) следует вычислить математическое ожидание тх ^) = М [ X ^)] и дисперсию Dx(0 = М[X2^)] векторной переменной X из линеаризованных методом статистической линеаризации уравнений (1) [9, 11]:

тх = б0

) = 2К, (тх, Dx, t)Dx + h\тх, Dx ^ (t), )х (to) = )х0,

где F0(mx,)х,t) - статистическая характеристика нелинейности ф - правой части уравнения (1), а Кй (тх,)х, t) - статистический коэффициент усиления.

Для гауссовой плотности вероятности /1(X, t) по формуле (8) при из со сторонами [—да, у] получим:

1

^ехр

(X — тх )2 2Л

dX.

Плотность вероятности /2(X, t) зададим усечённой гауссовой в следующем виде:

Г/20

^^ [0, X йи3,

где /20(X, t) - гауссова плотность вероятности с вектором математического ожидания тх ^) и дисперсией )'х ^), определяемыми из следующих уравнений на каждом г-м интервале

, ^ +Т3 ] [11]:

тх =ф0(тх,)х,t) — | XS(X,t)dX — q(t)mX,

т'х (ti) = тх ),

)X = 2КЕ (тХ, )Х, t))' + ^К, )Х )0Х (t) —

— | X2 3(X, t)dX — q(t))',

торой характеризуется уравнением (2), преобразованному к виду:

х = ах + x

где - белый шум интенсивности О иМ[^]= 0.

Область захвата из одномерная и определена так: —да < ф < у . Плотность вероятности распределения фазовой координаты х при гауссовом начальном условии также гауссова /1(х, t) с параметрами тх^), )х ^). Тогда вероятность Р^) в соответствии с (20) будет найдена по формуле:

У — тх

Pl(ti) = 1 — Ф

где Ф(...) - функция Лапласа, а тх(t) и )х(t) определяются из уравнений:

тх = атх , тх ) = тх0"

4 = 2аВх + О, )х

Плотность вероятности /2(х, t) в соответствии с допущением (21) представим в виде:

1

/2( X, t) = <

)'х (tг ) = )'х (tг ).

Интегрируя уравнения (22) и (23), с учётом формулы (21) и Р2(^-) = 1, определяем Р2(^- + + т3 ¡X(ti) е и3) по формуле (16).

После вычисления ) и Р2(^- +т3 [X (ti) е е и3) на каждом интервале определяем вероятность Р3 ) по формуле (17).

Пример

В качестве примера рассмотрим для наглядности одномерную систему, эволюция ко-

у12Пох

0, х > у.

ехр

(х — т'х )2

2Щ

, — да < х < у,

Для определения P2(t) в соответствии с (14) получаем уравнение:

Р2 = [—П2

где:

П2

2 дх

Для определения т'х ^) и )'х ^) запишем уравнения вида (22) и (23):

Рис. 2. График зависимости Р1

—да

0

—да

Р2 1

0.75

0.5

0.25

1

\ 1 м

\ I *

\ \ 4

\ * 1 ! * \

\ 1 ; 1

\ \

\ Ч 1 \

41 4 чI ч

0.1 0.2 0.3

Рис. 3. Графики зависимостей Р2

0.4

Рис. 4. Значения Р^) соединенные интерполяционной кривой

Р

2 п'

х '

т'х = атХ - xп2(У,t)т'х, т'х) = тх), (31) Р2

D'x = 2aD'x+G - (у - т'х)2 п (у, t) - РЦ

Р2

Ц'х 6) = Ц[ (tг ). Производную д^ (х, t) / дх на границе у вычисляем приближенно в соответствии с рекомендацией [8] следующим образом:

| [ f2( х, t )]х=у= А*^

дх ' Ах

где Ах - малое приращение координаты х.

Расчёты были проведены в среде Mathcad при следующих исходных данных: а = 3, О = 2, у = 0,5, Ах = 0,05, х3 = 0,1 с, % = 0,4 с.

На рис. 2 изображен график зависимости

Л(0.

На рис. 3 изображены графики Р2(0 на интервалах времени (^+1 - ti), 7 = 1...4.

На рис. 4 изображены значения Рз (^) = Р ^), г = 1-4, соединенные интерполяционной кривой.

Заключение

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

Приведённая выше методика позволяет приближённо решать ряд практических задач вероятностного анализа попадания БЛА в различные

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

слежения за БЛА по сравнению со временем пребывания БЛА внутри границ заданной области пространства (запретной зоне).

ЛИТЕРАТУРА

1. Беспилотные летательные аппараты [Электронный ресурс]. - Режим доступа: http://geoportal.by/index/ bespilotnye_letatelnye_apparaty/0-434/.

2. Семейство беспилотных авиационных комплексов видео-мониторинга местности «Бусел», «Бусел М», «Бусел М50» [Электронный ресурс]. - Режим доступа: http:// uavbusel.by/catalog.html/.

3. Roos, J. C. and Peddle, I. K. «Autonomous Take-off and Landing of an Unmanned Aerial Vehicle», R & D Journal of the South African Institution of Mechanical Engineering, volume 25, pp 16-28, 2009.

4. Моисеев, В. С. Прикладная теория управления беспилотными летательными аппаратами / В. С. Моисеев. - Казань: ГБУ РЦМКО, 2013. - 768 с.

5. Распопов, В. Я. Микросистемная авионика: учебное пособие / В. Я. Распопов. - Тула: «Гриф и К», 2010. -248 с.

6. Современные информационные технологии в задачах навигации и наведения беспилотных маневренных летательных аппаратов / под ред. М. Н. Красильщикова, Г. Г. Серебрякова. - М.: Физматлит, 2009. - 556 с.

7. Дмитриевский, А. А. Внешняя баллистика / А. А. Дмитриевский, Л. Н. Лысенко. - М.: Машиностроение, 2005. - 608 с.

8. Лобатый, А. А. Вероятностный анализ сенсорного автомата / И. Е. Казаков, А. А. Лобатый / Автоматика и телемеханика, 1986, N3, С. 74-79.

9. Пугачев, В. С. Теория стохастических систем / В. С. Пугачев, И. Н. Синицын. - М.: Логос, 2004. - 1000 с.

10. Robert J. Bamberger Jr., David P. Watson, David H. Scheidt, and Kevin L. Moore. Flight Demonstrations of Unmanned Aerial Vehicle Swarming Concepts / Johns Hopkins APL Technical Digest, Volume 27, Number 1 (2006). - 55 p.

11. Казаков, И. Е. Анализ систем случайной структуры / И. Е. Казаков, В. М. Артемьев, В. А. Бухалев. - М.: Наука, 1993. - 270 с.

REFERENCES

1. Unmanned aerial vehicles [Electronic resource]. - Access mode: http://geoportal.by/index/bespilotnye_letatelnye_ apparaty/0-434/.

2. The family of unmanned aerial systems for video monitoring of the area «Busel», «Busel M», «Busel M50» [Electronic resource]. - Access mode: http:// uavbusel.by/catalog.html/.

3. Roos, J. C. and Peddle, I. K. «Autonomous Take-off and Landing of an Unmanned Aerial Vehicle», R & D Journal of the South African Institution of Mechanical Engineering, volume 25, pp 16-28, 2009.

4. Moiseev, V. S. Applied control theory of unmanned aerial vehicles / V. S. Moses. - Kazan: GBU RCMKO, 2013. -768 p.

5. Raspopov, V. Y. Microsystem avionics: a training manual / V. Y. Raspopov. - Tula: «Grif and K», 2010. - 248 p.

6. Modern information technology in the problems of navigation and guidance of unmanned maneuverable aircraft / ed. M. N. Krasilshchikova G. G. Serebryakova. - M.: Fizmatlit, 2009. - 556 p.

7. Dmitrievsky, A. A. External ballistics / A. A. Dmitrievsky, L. N. Lysenko. - M.: Mechanical Engineering, 2005. - 608 p.

8. Lobaty, A. A. Probabilistic analysis of a sensory automaton / I. E. Kazakov, A. A. Lobaty / Automation and Remote Control, 1986, N3, S. 74-79.

9. Pugachev, V. S. Theory of stochastic systems / V. S. Pugachev, I. N. Sinitsyn. - M.: Logos, 2004. - 1000 p.

10. Robert J. Bamberger J.r., David P. Watson, David H. Scheidt, and Kevin L. Moore. Flight Demonstrations of Unmanned Aerial Vehicle Swarming Concepts / Johns Hopkins APL Technical Digest, Volume 27, Number 1 (2006). - 55 p.

11. Kazakov, I. E. Analysis of random structure systems / I. E. Kazakov, V. M. Artemyev, V. A. Bukhalev. - M.: Nauka, 1993. - 270 p.

Поступила После доработки Принята к печати

14.10.2019 29.11.2019 01.12.2019

LOBATY A. A., YATSYNA Y. F, STEPANOV V. Y, BUMAI A. Y.

PROBABILISTIC ANALYSIS OF GETTING INTO IN RESTRICTED AREA BY UNMANNED AERIAL VEHICAL

Belarusian National Technical University

The probabilistic analysis of crossing by an unmanned aerial vehicle (UAV) of the boundary of the no-fly area is solved. Condition of stating of the fact of the violation of the boundary of the restricted area is to stand of UAV within the area during

a specified time. The substantiation of the mathematical model to research through linearized vector stochastic equation is carried out. The problem is solved by applying the theory of Markov processes of random structure with absorption of realizations at the boundary of a given area. Particularity of the approach is the contemporaneously consideration of two probability densities of the distribution of phase coordinates that describe the boundary conditions. In this case, two equations systems are solved for probabilistic moments: taking into account the absorption of realizations and without taking into account the absorption. The probability of an object gets into specified area and do not leave one during the time that necessary to notice the unmanned aerial vehicle at the restricted area.

Keywords: unmanned aerial vehicle, boundary, structure, probability density, equations for moments.

Лобатый Александр Александрович, доктор технических наук, профессор. С 2000 года возглавляет кафедру «Информационные системы и технологии» Белорусского национального технического университета. Проводит исследования в области анализа и синтеза систем управления, в том числе - беспилотными летательными аппаратами. Автор и соавтор множества статей в научных журналах и конференциях, автор ряда книг и учебных пособий. Тел: + 375 (29) 346-82-56. E-mail: lobaty@bntu.by

Lobaty A. A., doctor of Science, Professor. From 2000 he heads of the department «Information Systems and Technologies» at the Belarusian National Technical University. Conducts research in the areas of analysis and synthesis of control systems including unmanned aerial vehicles. He is the author and co-author of many articles in scientific journals, conferences and books.

Яцына Юрий Францевич, директор Научно-производственного центра многофункциональных беспилотных комплексов Национальной академии наук Беларуси. Специалист в области разработки систем управления мобильных роботизированных систем. Руководит научно-исследовательской и производственной деятельностью по созданию беспилотных летательных аппаратов.

Yatsyna Y. F., Director of the Scientific-Production Center of the multifunctional unmanned systems of the National Academy of Sciences of Belarus. A specialist in the field of development of control systems of mobile robotic. Head of research and development work on the creation of unmanned aerial vehicles.

Степанов Владимир Юрьевич, исследователь в области технических наук, программист, преподаватель кафедры «Информационные системы и технологии» Белорусского национального технического университета. Научные интересы касаются: ИТ, программирования, робастного управления, обработки цифровых сигналов и параметрической оптимизации систем управления.

E-mail: vladimir.stepanov@bntu.by

Stepanov V. Y., Researcher in the field of technical sciences, Software Engineer, teacher of « Information Systems and Technologies» department of Belarusian National Technical University. His research interests include IT, programming, robust control, digital signal processing and parametric optimization of control systems.

Бумай Андрей Юрьевич, аспирант кафедры «Информационные системы и технологии» Белорусского национального технического университета. Проводит исследования в области анализа и синтеза стохастических систем управления применительно к беспилотным летательным аппаратам.

E-mail: andrei.bumai@gmail.com

Bumai A. Y., PhD student of «Information Systems and Technologies» department of Belarusian National Technical University. Conducts research in the areas of analysis and synthesis of stochastic control systems applying to unmanned aerial vehicles.