Марковская модель оценки риска катастроф на воздушном транспорте Текст научной статьи по специальности «Математика»

2005 НАУЧНЫЙ ВЕСТНИК МГТУ ГА №90(8)

серия Эксплуатация воздушного транспорта и ремонт авиационной техники.

Безопасность полетов

УДК 629.735

МАРКОВСКАЯ МОДЕЛЬ ОЦЕНКИ РИСКА КАТАСТРОФ НА ВОЗДУШНОМ ТРАНСПОРТЕ

В. Л. КУЗНЕЦОВ

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

Введение

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

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

При описании перемещения ВС в пространстве используются различные пары векторных функций, характеризующих: реальные положение и скорость ВС как материальной точки (r(t ),v(t)), плановые параметры полета (r0(t),v0(t)), а также координаты и скорость объекта,

зафиксированные средствами наземных служб или на борту ВС (f ,vf).

Центральным моментом в задаче оценки риска катастроф является вопрос о виде и способе задания функции распределения плотности вероятности состояния ВС в фазовом пространстве

- F (t, r ,v). Поэтому мы будем различать две постановки задачи.

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

Fi (A r,v ) = /г 4 (t ),V0i (t)) , (1)

где i - индекс, приписываемый ВС. В модели разумно ограничиться случаем двух ВС, т. е. i = 1,2. Отметим, что в этой постановке задачи зависимость F от t содержится в плановых параметрах полета.

Во второй более сложной постановке задачи распределения случайных отклонений считаются известными лишь в некоторые дискретные моменты времени - {tk},к = 1,2,..., а необходимая для вычисления риска катастрофы функция распределения

Fi (t, r,V ) = f (t, r,Vfi (tk ),Vf г (tk )) t > tk

находится как решение задачи Коши для кинетического уравнения

Э F - Э F - dFi n

~dt~ + V ~ЭТ~ + 1; ~ЭА = 0, tk < t < tk+1 . (2)

Здесь в качестве начального условия берется функция

Рг (гк,Г,У) = /■ (*,г,у\ цг {гк), ^ {гк)) _ =Фг {г,у\ цг {гк), ^ (гк)),

г к

вид которой определяется ошибками измерения состояния ВС и полагается известным .

В простейшем случае отсутствия нерегулярных внешних воздействий и управления ВС

( V _ 0 ) решение (2) получается просто: ^ (г, ГУ) _ Фг Г + V' г, V | г}г (гк ), (гк )).

В наиболее интересном случае, при учете управления, компенсирующего отклонения ВС от

плановых параметров полета (V Ф 0), кинетическое уравнение (2) следует дополнить динамической моделью ВС, учитывающей конкретные особенности управления (режим

автопилота, ручное управление и т.д.) и случайные внешние воздействия (например, поле ветровых потоков).

Инвариантным ядром этих двух постановок является задача об оценке риска катастрофы при заданном для ВС распределении плотности вероятности в 6-мерном фазовом пространстве (г ,у ) как функции времени. По существу именно такую задачу мы будем рассматривать в этой работе. Однако, для наглядности изложения, ее удобно представить в виде первой постановки, а именно: в некоторой области трехмерного пространства движутся два ВС, состояния которых описываются зависящими от времени функциями распределения (1) в фазовом пространстве (г, V). В начальный момент времени, при г _ 0, достоверно известно, что катастрофы не было. Необходимо вычислить вероятность катастрофы на временном интервале (0, г), т.е. найти функцию Рст _ Рст (г).

Очевидно, что Рст (г) - неубывающая функция времени, качественный вид которой представлен на рис.1.

на временном интервале (0, г )

Как будет видно из дальнейшего, вероятность катастрофы на временном интервале (г1, г2), при условии, что до момента времени г1 столкновения воздушных судов не было, равно ЛР _ Рст (г2) _ Рст (г1) _ Р2 - Р . В задаче оценивания уровня безопасности риском катастроф называют величину, равную предполагаемому числу столкновений ВС за единицу летного времени (1час). Если отнести это определение к задаче возможного столкновения двух ВС, то интервал Лг _ г2 - г1 _ 1 час - достаточно большая величина и, как видно из рисунка, на этом

интервале Рст (t) претерпевает как быстрые изменения (период максимального сближения ВС), так и имеет области, где Рст (t)» const. Деление на At = 1 час фактически вуалирует опасную область (t[, 12), уменьшая величину среднего риска. Поэтому целесообразно ввести понятие локального риска, определяемого по формуле

R(t) = Aim Рт(t + AA)~Рт(<) = С) • (3)

At®° At dt

и оценивать риск столкновения на интервале (°, t) как

R = max R(t), (4)

tG (°,t)

т.е рассматривать в качестве меры риска максимальный локальный риск, оцененный по формуле (3).

Из сказанного видно, что основной проблемой в задаче оценивания уровня безопасности является вычисление Рст (t), в основе расчета которой лежит аналитическое задание

распределения отклонений состояния ВС от номинальных (плановых). Наиболее часто цитируемыми в такой постановке являются модели Рейха [1] и Хсю [2]. В настоящей статье предлагается модель, объединяющая и корректирующая подходы этих авторов.

1. Описание столкновения ВС с помощью марковского процесса

Введем случайную функцию X = X(w, t), t g [°,Г], wg W (W - пространство элементарных событий), равную нулю, если до момента времени t столкновения не произошло и единице в противном случае. w- случайное событие, при фиксировании которого (w= w) t) становится неслучайной функцией t , называемой реализацией случайного процесса. Пример такой реализации показан на рис.2. Если зафиксировать другое элементарное событие , то скачек графика, представленного на рис.2, будет иметь место уже при другом значении t, или его вообще не будет - катастрофы не произойдет. Различные а соответствуют различным возможным траекториям ВС, реализуемым при фиксированных плановых

параметрах полета - (r0i (t),v° (t)) i = 1,2 .

Рис.2. Вид одной из реализаций случайного процесса X = і ). і - момент катастрофы

Очевидно, что для процесса X = %(®, і) выполняется марковское условие, т.е. для любых

і1 , і2, і3 є [0, Т] и удовлетворяющих условию 0 < і1 < і2 < і3 выполняется соотношение

Р^З )| 2 ), x(t1 )] = )| 2 )]

(5)

Действительно, состояние ВС в момент времени ¿3 предопределено, если %(2) = 1 (к моменту ¿2 катастрофа произошла). И при этом неважно, что было с ВС в момент времени ¿1. Если %( 2) = 0, то XI) тоже должно быть равно нулю (до ¿2 катастрофы не было), и поэтому задание Х(А) - избыточная информация.

Пусть далее Р1^) - вероятность того, что ) = 0, Р2(^)- вероятность противоположного

случая. Эволюция этих вероятностей описывается уравнениями Колмогорова [3], которые в нашем случае имеют вид:

— = -12 ' Р1

Жі ,2 1

жр2

Жі

(6)

= Л,2 • р

Начальные условия для нашей задачи задаются соотношениями:

Р(0) = 1, Р? (0) = 0 .

В (6) Л12 - плотность вероятности перехода системы из состояния X = 0 в состояние

£ =1

Л,2(і )

ІІШ

Аі ——0

р£(і + А) = 1 Х(і) = о)

Аі

(7)

Здесь р(§(4 + Ьл) = 11 X) = 0) - условная вероятность того, что к моменту времени ( + Ь)

столкновение произойдет, при условии, что в момент времени I катастрофы еще не было. Решение (6) для Р2 (I) имеет вид

р2(і) = 1 - ехрі |1,2(Т)Жт!

(8)

хорошо известный в литературе [1,2,4,5] и его можно было бы не приводить, если бы не трактовка соотношения (7) как плотность вероятности перехода между состояниями.

В широко цитируемой работе Хсю[2] и более поздних статьях [4,5] эта величина определяется формулой (см., например, (3) в [2])

1(і) =

плотность вероятности перекрытия ВС в мом. і объем области столкновения

линейный размер области перекрытия

относительная скорость ВС]

(9)

Некорректность этой формулы заключается в том, что она основывается на использовании безусловной, а не, как следует из (7), условной вероятности катастрофы. Более последовательный подход реализуется в модели Рейха [1], использующей формулу Райса для потока векторного поля через поверхность. Однако, в отличие от модели Хсю, где реализуется переход в систему отсчета, связанную с плановым положением ВС1 («регулярная» система отсчета), в работе Рейха новая система связывается с реальным положением ВС1 и в этом смысле ее положение относительно земной - случайно. Вероятно, сложности пересчета

о

вероятностных распределений в такие системы отсчета и стимулировало появление работы Хсю [2].

Отметим здесь также, что отсутствие предельного перехода по At мало сказывается на величине результата, но позволяет перейти от объемных интегралов, которые фактически присутствуют в числителе (9) у Хсю [2], к поверхностным, описывающим потоки в модели Рейха [1]. Подробнее эти вопросы рассматриваются во втором разделе нашей статьи и в приложении.

Вычисление плотности вероятности перехода между состояниями

Рассмотрим пару ВС плановые перемещения которых описываются вектор-функциями r01(t) и r02(t) в земной системе координат. Перейдем в систему, связанную с первым ВС (ВС1).

Отметим, что поскольку истинное положение ВС1 - случайно, то отдельно возникает вопрос о вычислении ошибок положения ВС2 в этой новой системе отсчета. Этот вопрос мы рассмотрим в разделе 3, а сейчас будем полагать, что распределение ошибок состояния ВС2 известно и задается плотностью распределения w(r',v'2|r01(t), r02(t)). Здесь штрихованные величины соответствуют новой системе координат. Плановые характеристики полета, исполняющие роль зависящих от времени параметров распределения, в новую систему не пересчитываются. V2 -реальная скорость ВС2 в системе, связанной с ВС1 (по модулю совпадает с относительной скоростью движения ВС).

Пусть G - критическая область, т.е. область вокруг ВС1, при попадании центра ВС2 в которую неизбежно произойдет катастрофа. Далее мы будем полагать G сферой радиуса R, хотя полученные результаты справедливы для любой выпуклой области, в частности для цилиндра.

Событие A, заключающееся в том, что V2 е [W,W + AW] и к моменту времени (t + At) катастрофа произойдет при условии, что до момента времени t столкновения не было, может иметь место в том, и только в том случае, если в момент времени t ВС2 было вблизи области G (область dG на рис.3) но не внутри нее, как предполагает Хсю [2]. Действительно, находясь

в dG в момент времени t и обладая скоростью W, ВС2 за интервал времени At «успевает» попасть в область G -катастрофа происходит.

Если предположить, что w(r2, V2|r01(t), r02(t))» const внутри областей G и dG, а предел в

(7) можно убрать, ограничившись некоторым малым, но фиксированным значением At, то специальный выбор интервала времени позволяет сделать равными интегралы по обеим указанным на рис.3 областям. Именно этим обстоятельством и воспользовался Хсю, выбрав в качестве At среднее время «пролета» ВС2 области G .

В соответствии со строгой формулой (7) At необходимо устремить к нулю. Тогда для малых At dG есть множество точек, удаленных от границы G на расстояние не большее, чем

|(V2, n)| - At, где n - внутренняя нормаль к границе dG области G в точке ее пересечения

реализацией случайной траектории ВС2. Вероятность события A может быть записана в виде:

Р(A) = f W(r^V2 |?01 (t),Г02 (t)) - AW »

dG f . (10)

At f (V2,n )- q((V2,n ))w(r2,V2 Ir01 (t),r02 (t)) ds - AW

d G

Рис.3. Графическая иллюстрация соотношения области столкновения по Хсю (О) и области грядущего столкновения (dG), определяющей переход в состояние катастрофы за время Аі

Проинтегрировав (10) по спектру возможных относительных скоростей находим:

р(х(і + Аі) = ^ Х(і) = 0) = Аі|^ • І(т^Л)• в(К,ЛК1(і)/02(і))(11)

до

Отсюда для 12 (і) имеем

) _ | ЖТ2 • | (у,2,Л )^((У,2,Л |^01(і),^02(і)) . (12)

до

В (10)-(12) в((у'2, Л)) - в-функция Хевисайда, учитывающая тот факт, что ВС может только влететь в критическую область, но не вылететь из нее. Т.е. в поверхностном интеграле (12) следует учитывать только «входящий» в О поток плотности вероятности. Это соответствует символу «+» в модели Рейха.

Для получения решения задачи соотношения (8),(12) следует «доукомплектовать» явным видом функции П’(г2', Т'2 |г01(і), г02(і)). Перейдем к обсуждению этого вопроса.

2. Вычисление плотности вероятности состояния ВС2 в собственной системе отсчета ВС1

В этом разделе будем полагать, что известны распределения плотности вероятности состояний обоих ВС в шестимерном фазовом пространстве, связанным с землей, т.е. известны функции распределения /1 (г,у| г01(і)) и /2(г,у| г02(і)). Необходимо получить вероятностное

описание ВС2 в системе отсчета, связанной с ВС1, положение и скорость которого случайны.

Изменение состояния (6-вектора (г (і),у(і))) ВС при переходе из неподвижной (земной) системы отсчета в движущуюся описывается преобразованием Галилея:

г'(і') = г (і) - V • і, і ' = і (13)

исследующему из них закону сложения скоростей

у' (і) = у(і) - V (14)

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

При переходе в собственную систему координат ВС1 соотношения (13),(14) приводят к следующему выражению для вектора состояния ВС2 в фазовом пространстве

(гу) = (г2- /■, у2 - V!) (15)

Из (15) следует, что вектор состояния (г', V ) представим в виде разности двух случайных 6мерных векторов (г2, У2 ) и (г , V ) .

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

н(г', 0 X 4 0)) = | ^ ^ ¡1 (г1 д |г01 (г ))/^2 (Г+г1, у' +у |г02 ^)) (16)

Формулы (8),(12), (16) дают решение поставленной задачи, т.е позволяют определить вероятность катастрофы двух ВС, перемещающихся по плановым траекториям г01 (^) и г02 (^) .

В Приложении приводится наше видение вывода основных соотношений работы [2], которое, совместно с замечаниями к формуле (9), позволяет рассматривать модель Хсю, как огрубленную модель Рейха. У Хсю вычисление поверхностного интеграла второго рода для плотности потока вероятности (12) заменяется объемным интегралом и затем по теореме о среднем - простым произведением плотности вероятности перекрытия на объем области столкновения (9).

Заключение

В работе предлагается марковская модель расчета риска катастроф, развивающая модель Рейха [2]. Показано, что модель Хсю [2] можно рассматривать как огрубленную, но удобную для расчетов модель Рейха.

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

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

Приложение: вывод выражения для плотности вероятности перекрытия в модели Хсю

Рассмотрим реальное взаимное расположение двух ВС в системе координат, связанной с плановым (номинальным) положением ВС 1 (рис. 4).

Рис. 4 - Схема взаимного расположения ВС в модели Хсю

На рис. 4 введены следующие обозначения: r1 ■ радиус-вектор, описывающий реальное положение BC1, S0 (t) - радиус-вектор планового положения BC 2, r2 - вектор,

описывающий отклонение в положении BC 2 от планового, р - вектор смещения геометрического центра BC 2 относительно от центра BC 1, GR (r1 ) - критическая область перекрытия с центром в точке реального расположения BC 1 .

Из рис. видно, что

Г2 + S o(t) = ri +Р (П. 1)

Если |р| < R - т.е. конец вектора r2 принадлежит области GR, то имеет место столкновение.

Полагая реальные положения BC независимыми и распределенными в соответствии с плотностями распределения f1(r1), f2(r2), можно найти вероятность столкновения Рст(t) в момент времени t с помощью следующих рассуждений: вероятность того, что центр BC 1 будет находиться в окрестности точки r - ( f1(r,)^r1), а центр BC 2 попадет в критическую область ( J f2 (Г2 )dr2 ) в силу независимости положений BC определяется произведением

Gr (r-So(t ))

f1(r1)dr1 • J f2(r2)dr2 . (П. 2)

Gr (r1 -So(t ))

Интегрируя (П. 2) по возможным реальным положениям BC 1 , находим

Pcm (t ) = jÆJlifi) j f2 )d2 (Д 3)

2

ок (Г -£,(г))

Фигурирующий в (П. 3) интеграл является повторным, т.к. область интегрирования внутреннего интеграла зависит от г1 .

Делая последовательно две замены переменных во внутреннем интеграле (П.3) г2 = Г + р - £0 (г), £0 - р = £ , приходим к выражению:

Рст (0 = | ЯМА) | Л(А1 - £0(0 +Р)йр = | <ЗГ1/1(Г1) | /2 (Гг- А ^ (П. 4)

Оя (°) Ок (§0(1))

Теперь в повторном интеграле в (П.4) можно изменить порядок интегрирования. В результате получаем:

Р, (г) = | Я - С (А) , (П. 5)

Ря (*А0 (г))

где С(А) = | ^ (А )М2 (А - А) - введенная Хсю плотность вероятности перекрытия ВС.

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

ЛИТЕРАТУРА

1. Reich, P.G. Analysis of long-range air traffic systems - separation standards. Part I, II, III // The Journal of Institute of Navigation, 1966, vol.19, p.88,169,331

2. Hsu, D.A. The evaluation of aircraft collision probabilities at intersecting air routes // The Journal of Navigation, 1981, vol.34, N1, p.78

3. Розанов Ю.А. Теория вероятностей, случайные процессы и математическая статистика.

М., Наука, 1985

4. Lezaud P., Mehadhebi K A synthesis of current collision risk models, SASP/WHL/5/WP30, Tokyo, Japan, 17-28 May 2004

5. Грибков И.М., Спрысков В.Б., Щербаков Л.К. Модель оценки риска катастроф ВС при движении по пересекающимся воздушным трассам на одной высоте.// Научный вестник МГТУГА. Настоящий выпуск.

6. Королюк В.С., Портенко Н.И., Скороход А.В., Турбин А.В. Справочник по теории вероятностей и математической статистике. М., Наука, 1985

MARKOV COLLISION RISK MODEL ON AIR TRANSPORT

Kouznetsov V.L.

Markov collision risk model on air transport as a development of Reich model is suppOsed. The constraints whereby Reich model reduced to Hsu model are discussed.

Сведения об авторе

Кузнецов Валерий Леонидович, 1949 г.р., окончил МГУ им. Ломоносова (1972),

доктор технических наук, заведующий кафедрой прикладной математики МГТУ ГА, автор более 80 научных работ, область научных интересов - методы математического моделирования в задачах распространения излучения в пространственно неоднородных случайных и периодических средах.