«К. К. ВАСИЛЬЕВ, М. Н. СЛУЖИВЫЙ МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ СИСТЕМ СВЯЗИ Учебное пособие по дисциплине Математическое моделирование каналов и систем телекоммуникаций для студентов специальностей 21040665 Сети связи и ...»

Модели СП в виде стохастических ДУ обладают рядом положительных свойств. Прежде всего, в отличие от других методов описания СП (например, с помощью ПРВ или моментных функций), они непосредственно определяют способ генерации его реализаций, для осуществления которого можно вполне успешно использовать аналоговые вычислительные устройства или соответствующие программы для ЭВМ. Использование стохастических ДУ в форме уравнений состояния позволяют синтезировать алгоритмы оптимальной обработки и генерации для широкого класса сигналов и помех [16, 32, 44].

Для большинства реальных СП используется модель в виде ДУ, линейных относительно входных воздействий типа белого шума где f ( x,t ) - векторная и v ( x,t ) - матричная функции векторного аргумента [16, 44]; (t ) - белый гауссовский шум с нулевым математическим ожиданием и заданной корреляционной функцией (КФ). Задавая различные функции f ( x,t ) и v ( x,t ) в уравнении (3.2), можно получить СП с различными статистическими характеристиками. Так, если f ( x,t ) – линейная функция, а v ( x,t ) не зависит от x (t ), то уравнение (3.2) будет определять гауссовский СП. В работах [16, 38] приведены связи между функциями f ( x,t ), v ( x,t ) и соответствующими статистическими характеристиками (ПРВ, КФ).

Следует отметить, что задача синтеза марковской модели СП, т. е. определение функций f ( x,t ) и v ( x,t ) по заданным статистическим характеристикам, не всегда имеет одно решение. Здесь возникает вопрос о единственности решения. Если заданы статистические характеристики моделируемого процесса в форме функций распределения, то можно указать несколько уравнений различных типов, порождающих процессы с такими характеристиками. В случае нелинейных моделей это является неизбежным следствием ограниченности исходных сведений о процессе. Если единственное решение задачи синтеза отсутствует, то при выборе уравнения, которое будет заложено в основу функционирования какого-либо устройства, можно учитывать соображения сложности, экономичности и так далее. На практике число эквивалентных, с точки зрения статистических характеристик, уравнений обычно не велико (не более трех, четырех), и из них всегда можно выбрать единственное [32].

В последнее время, в связи с широким использованием цифровых устройств в радиотехнических системах для статистического описания сигналов и помех используют случайные последовательности, которые могут быть заданы разностными стохастическими уравнениями где xi - состояние системы в момент времени ti ; ( xi 1,i ) - векторная функция;

v ( xi 1,i 1) - матричная функция; i – последовательность СВ [32, 33]. Если в качестве возбуждающей последовательности i использовать дискретный белый шум, то уравнение (3.3) будет определять марковскую последовательность.

При использовании стохастических уравнений (3.1), (3.3) в качестве моделей случайных сигналов и помех возникают задачи анализа и синтеза [38].

Первая из них заключается в определении статистических характеристик (в первую очередь ПРВ и КФ) СП или случайной последовательности по заданным функциям f, v и. Более сложной задачей синтеза является определение неизвестных функций f, v и по заданным статистическим характеристикам.

Для непрерывных СП обе эти задачи решены для достаточно большого числа частных случаев [8, 16, 23]. При этом анализ выполняется на основе решения уравнения Фоккера-Планка-Колмогорова [44].

Для случайных последовательностей эти задачи до конца решены лишь для случая линейной функции ( xi 1,i ) и независящей от xi 1 функции v [5].

Кроме того, в работах [32, 38, 44] рассмотрено решение задачи синтеза и для нелинейной функции. Однако при этом накладываются ограничения на вид ПРВ. Отсутствие конструктивных решений задачи синтеза, применительно к разностному уравнению (3.3), не позволяет получить ММ случайных последовательностей с заданными статистическими свойствами.

3.3. Нелинейные модели марковских случайных процессов В данном разделе рассматривается задача построения стохастических моделей негауссовских СП с заданной ПРВ. Возможность достаточно простого определения нелинейных функций f ( x,t ) и v ( x,t ) для стохастических уравнений (3.2) обусловлена существованием обыкновенного ДУ, связывающего эти функции с ПРВ w ( x ) при t. При этом задача нахождения нелинейных функций f ( x,t ) и v ( x,t ) наиболее проста если ПРВ w ( x ) входит в класс распределения Пирсона [44]. Вместе с тем, для стохастических разностных уравнений (3.3) лишь для линейных систем может быть точно определено безусловное распределение.

В работах [16, 38] подробно исследованы стохастические ДУ, определяющие различные негауссовские процессы. Запишем векторное стохастическое ДУ, описывающее непрерывную динамическую систему где N ( t ) - диагональная матрица с неотрицательными элементами, ( t ) - стандартный винеровский процесс с M {d ( t ) d T ( t )} = E dt ; E - единичная диагональная матрица. При этом стационарная ПРВ будет:

(3.5) (3.6) 3) логарифмически-нормальной (3.7) Стохастические ДУ (3.4)-(3.7) исследованы наиболее полно, дают адекватное описание СП во многих радиотехнических системах и позволяют получить замкнутые выражения для алгоритмов оценивания сигналов. Однако системы, базирующиеся непосредственно на таких уравнениях, как правило, не могут быть реализованы и исследованы из-за отсутствия необходимых стабильных нелинейных элементов аналоговой техники. В связи с этим на практике осуществляется переход от непрерывных к дискретным во времени системам на основе различных разностных схем [32, 33]. Для обеспечения соответствия разностной схемы и уравнений вида (3.3), (3.4) от ЭВМ требуется высокое быстродействие и большая точность представления чисел.

Стохастические разностные уравнения (3.3) свободны от этого недостатка, так как представляют последовательность чисел, для формирования которой можно воспользоваться методами цифровой техники. Поэтому в задачах с нелинейными моделями наблюдений, но при гауссовских компонентах, системы с дискретным временем имеют очевидные преимущества при технической реализации. Вместе с тем, при отличии ПРВ случайных последовательностей (3.3) от гауссовских возникает задача определения нелинейных функций ( xi 1,i ) и v ( xi 1,i 1) по заданным распределениям.

В ряде задач желательно построить (синтезировать) нелинейное уравнение (3.3) на основе априорных данных о виде безусловного одномерного распределения. К сожалению, такого же простого пути, какой существует для построения моделей (3.5)-(3.7), здесь нет. В работе [15] синтез марковской модели приводится не на основе априорных статистических характеристик СП, а на экспериментальных данных. При этом решение получено лишь для ограниченного класса линейных функций ( xi 1,i ), v ( xi 1,i 1). Полученные модели марковских СП используются для решения задач прогноза.

Рассмотрим асимптотический подход для решения поставленной задачи, основанный на предположении о существовании предельного стохастического ДУ Ито [16, 38] для заданного разностного стохастического уравнения (3.3), записанного в виде Будем связывать значения xi 1 и xi c моментами времени ti1 = t и ti = t + t. При условии существования предельных соотношений для детерминированных функций перепишем (3.9) в форме, позволяющей найти характеристики приращений xi xi 1 :

нормальному закону с нулевым математическим ожиданием и единичной дисперсией.

3.3.1. Моделирование случайных последовательностей с ПРВ семейства распределений Пирсона Каждому конкретному распределению из семейства распределений Пирсона, которое описывается с помощью ДУ соответствует свой набор коэффициентов a0, a1, b0, b1, b2. Значит, задавая конкретное распределение, мы тем самым задаем и коэффициенты.

В качестве примера рассмотрим гамма-распределение, для которого ПРВ запишется в виде [44] где и - параметры распределения, Г( • ) - гамма-функция [7]. Вычисляя производную в (3.11), получим Сравнивая (3.11) и (3.10), найдем a0 =, a1 = 1, b1 =, b0 = b2 = 0.

Рекуррентное разностное уравнение запишется в следующем виде [44] Учитывая слагаемые более высокого порядка малости, вместо (3.13) получим Для конкретных значений и величина Vx может быть выражена через них, т.е. Vx = ( + 1) и уравнение (3.13) будет содержать только один параметр t.

полученной с помощью алгоритма (3.14) при t =0. Рис. 3.3. График теоретической ПРВ (сплошная линия) и ПРВ, полученной Для экспериментальной проверки степени совпадения ПРВ последовательности (3.14) с гамма-распределением было проведено статистическое моделирование синтезированного алгоритма. На рис. 3.2 и 3. представлены теоретическое распределение и распределения, полученные в результате моделирования. Весь диапазон изменения xk был разбит на одинаковых интервалов. Экспериментальные распределения получены на основе независимой выборки объемом 100000. Проведена проверка гипотез о принадлежности экспериментальных распределений гамма-распределению с помощью критерия. Для заданного числа интервалов (100) и выбранного уровня значимости ( = 0.005 ) критическое значение =140. При значениях t 0.003 некоторые реализации (3.14) приводили к неустойчивости. Для устранения этого явления, по всей видимости, необходимо использовать приближения более высокого порядка.

Динамическая система, описываемая разностным уравнением (3.14), является нелинейной. Поэтому провести более подробный анализ устойчивости этих алгоритмов не представляется возможным.

3.3.2. Моделирование случайных последовательностей Распределение Накагами часто используется для описания амплитудных замираний радиосигнала, прошедшего через турбулентную среду. ПРВ Накагами записывается следующим образом:

Нелинейное стохастическое ДУ, описывающее СП с ПРВ Накагами имеет вид В работе [44] получен рекуррентный для моделирования случайных последовательностей с ПРВ Накагами (при m = 0.5 ) Для проверки степени совпадения ПРВ последовательности (3.15) с распределением Накагами было проведено статистическое моделирование синтезированного алгоритма [16, 44]. При моделировании случайной последовательности с помощью алгоритма (3.15) удовлетворительные результаты были получены только для следующих значений коэффициентов:

m = 0.5, N = 1, a = 1. Попытка моделирования при t >0.5 оказалась неудачной, динамическая система при любом выборе коэффициентов a, N оказывалась неустойчивой. Полученные распределения значительно отличались от теоретической ПРВ. При получении случайных последовательностей с ПРВ Накагами с помощью алгоритма (3.15) можно изменить способ определения коэффициентов a и N.

Рис. 3.4. График теоретической ПРВ (сплошная линия) и ПРВ, полученной Рис. 3.5. График теоретической ПРВ (сплошная линия) и ПРВ, полученной На рис. 3.4 и 3.5 представлены теоретическое распределение и распределения, полученные в результате моделирования. Весь диапазон изменения xk был разбит на 100 одинаковых интервалов. Экспериментальные распределения получены на основе независимой выборки объемом 100000.

Проведена проверка гипотез о принадлежности экспериментальных распределений семейству распределений Накагами с помощью критерия.

Для заданного числа интервалов (100) и выбранного уровня значимости ( = 0.005 ) критическое значение =140. При значениях t 0. некоторые реализации (3.15) приводили к неустойчивости. Для устранения этого явления, по всей видимости, необходимо использовать приближения более высокого порядка.

3.4. Модели случайных процессов в виде временных рядов 3.4.1. Авторегрессионные модели случайных последовательностей Пусть последовательность СВ удовлетворяет стохастическому уравнению с начальным условием x0 = 0 ( 0 ), где и 0 – некоторые функции;

= ( 0, 1, 2,...) – заданная последовательность независимых СВ, называемая порождающей или возмущающей последовательностью.

Уравнение (3.16) представляет собой простейшую авторегрессионную модель СП X = { x0, x1, x2,...} [20].

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

Примером марковской СП может быть процесс авторегрессии 1-го порядка, последовательности 1,2,...,i,... независимых гауссовских СВ { i } по следующему правилу:

где < 1; M { i } = 0; M { i2 } =. Каждое очередное значение x i содержит часть предыдущего xi 1 и добавку в виде независимой СВ i [2, 5].

На рис. 3.6 представлены типичные графики реализаций такого процесса при различных значениях параметра, входящего в модель (3.17). Во всех случаях параметр, влияющий только на масштаб по оси ординат, выбран равным единице. Из этих рисунков видно, что при, близких к единице, процесс становится более гладким; при малых, напротив, значения процесса слабо зависимы между собой; при отрицательных корреляция между соседними значениями процесса отрицательна, поэтому он часто меняет знак.

При выборе начального значения x1, обеспечивающего стационарность и постоянство дисперсии M {xi2 }= x, параметр равен коэффициенту корреляции между любыми двумя соседними значениями СП. Действительно, умножая левую и правую часть (3.17) на xi 1 и находя математическое ожидание, получим M { xi xi 1} = 2 или = M { xi xi 1} 2.

Повторяя аналогичные операции после подстановки в уравнение (3.17) xi 1 = xi 2 + i 1, xi 2 = xi 3 + i2..., можно записать следующую формулу для КФ:

Таким образом, СП (3.1) имеет экспоненциальную КФ. В то же время СП (3.17) является марковской, поскольку любые вероятностные характеристики значения x i полностью определяются только предшествующим значением СП x i 1. При заданном x i 1 формула (3.17) позволяет найти все характеристики x i без учета предыстории, т. е. значений x1, x2,..., xi 2 СП. Так, условная ПРВ i = xi xi 1 и правил нахождения ПРВ функций СВ. Легко записать выражение и для совместного распределения произвольного числа n членов рассмотренной марковской СП:

Поскольку вид всех ПРВ перехода (3.18) не зависит от номера члена СП, Для стационарности необходимо выбрать СВ x1 таким образом, чтобы все безусловные ПРВ w ( x i ), i = 1, 2,..., n были одинаковыми. Проведенный анализ последовательности имеют нулевое среднее и дисперсию 2 = (1 2 ).

Кроме того, CП {x i } гауссовская, так как получена в результате линейного преобразования (1.42) гауссовских СВ {i }. Таким образом, ПРВ всех значений стационарной последовательности (3.17) будут иметь следующий вид:

x1 формируется как нормальная СВ с нулевым средним и дисперсией x, а последующие члены последовательности образуются в соответствии с рекуррентным соотношением (3.17).

Уравнения вида (3.17), которые часто называются уравнениями авторегрессии или стохастическими разностными уравнениями, представляют весьма узкий класс гауссовских марковских СП с экспоненциальной КФ.

Вместе с тем имеются различные возможности для существенного расширения этого класса [2, 5, 20]. Одной из них является описание СП с помощью авторегрессионных (АР) уравнений более высокого порядка:

где m порядок авторегрессии. С помощью подбора коэффициентов с разнообразными корреляционными свойствами [2]. Действительно, умножая (3.19) на xik и находя математические ожидания, получим после деления на x2 = M {x i2 }, i = 1, 2,..., n, следующее соотношение для значений КФ:

Общее решение этого разностного уравнения в стационарном случае представляется суммой экспонент [33]:

где v = ln zv ; zv, v = 1, 2,..., m, - корни характеристического уравнения z m 1 z m1 2 z m2... m = 0. Требование стационарности СП (3.20) выполняется, если v < 0, т.е. когда все корни zv, v = 1,2,..., m характеристического уравнения лежат внутри единичного круга на комплексной плоскости.

Подставляя в (3.21) значения k = 1, 2,..., m получим известную систему уравнений Юла-Уокера [2, 5]:

Решение этой системы позволяет найти коэффициенты 1, 2,..., m уравнения авторегрессии (3.20) по заданным или оцененным на основе эксперимента значениям Rx (1), Rx ( 2 ),..., Rx ( m ) КФ СП.

В качестве примера рассмотрим процесс авторегрессии второго порядка:

xi = 1 xi 1 + 2 xi 2 + i, i = 2, 3,..., n. Для стационарности процесса необходимо, чтобы корни характеристического уравнения z 2 1 z 2 = лежали внутри единичного круга, т. е. чтобы параметры 1 и 2 находились в треугольной области, показанной на рис. 3.7 [2, 5].

Рис. 3.7. Область значений коэффициентов 1 и 2 стационарного СП Значения КФ стационарной СП связаны между собой рекуррентным соотношением Rx ( k ) = 1Rx ( k 1) + 2 Rx ( k 2 ), k > 0, с начальными условиями Rx ( 0) = 1 и Rx (1) = 1 (1 2 ). Из этого соотношения следует, что Уокера 1 + 2 Rx (1) = Rx (1), 1Rx (1) + 2 = Rx ( 2 ) позволяет определить коэффициенты 1 и 2 уравнения авторегрессии по заданным или измеренным значениям R x (1) и R x (2 ) КФ.

Вид КФ определяется областью треугольника допустимых значений характеристического уравнения действительны и КФ представляет сумму двух затухающих экспонент. При 1 > 0, 2 > 0 (область I на рис. 3.7) корни имеют разные знаки:

осциллирующее слагаемое A2 z2 = A2 ( 1) z в области I коэффициент A1 > A2 и КФ R x (k ) не изменяет знака. Во второй области, показанной на рис. 3.7, оба корня положительны и КФ монотонно убывает. На одной границе области II ( 2 = 0 ) авторегрессия имеет первый порядок и Rx ( k ) = 1.

На другой границе, где 1 + 42 = 0 характеристическое уравнение имеет кратный корень z1 = z2 = 1 2. В этом случае выражение для КФ запишется области рис. 3.7 корни характеристического уравнения комплексные и КФ определяется по следующей формуле: Rx ( k ) = A0 d sin ( 0 k + ), где этом графики КФ имеют вид синусоиды с экспоненциальным уменьшением амплитуды.

Для иллюстрации рассмотренных ситуаций на рис. 3.8 и рис. 3. представлены зависимости КФ R x (k ) при различных значениях параметров и 2 АР уравнения. При построении зависимостей КФ на рис. 3. коэффициенты ( 1, 2 ) подбирались из различных областей треугольнике допустимых значений (рис. 3.7), но с учетом дополнительного условия R x (1) = M {x i x i 1 } = 0,9. Для всех КФ, представленных на рис. 3.9, таким дополнительным условием является один и тот же интервал корреляции k 0 = на уровне 0,5, т.е. R x (k = 7 ) = 0,5 [5].

3.4.2. Модели авторегрессии-скользящего среднего Стохастический линейный процесс можно представить как выходной сигнал линейного фильтра, на вход которого поступает белый шум i (рис. 3.10) передаточной функцией фильтра [2].

Последовательность 1, 2,..., образованная весами, теоретически может быть конечной или бесконечной. Если эта последовательность (конечная или бесконечная) сходящаяся, фильтр называется устойчивым, а процесс xi будет стационарным.

Рис. 3.10. Представление временного ряда с помощью линейного фильтра Модель авторегрессии (3.20) выражает отсчет xi процесса в виде конечной взвешенной суммы n предыдущих отсчетов процесса xi 1, xi 2,K, xi n плюс случайный отсчет i. Другой тип моделей, имеющий большое значение в описании СП, – это так называемый процесс скользящего среднего. Пусть xi линейно зависит от конечного числа m предыдущих отсчетов :

Такой процесс называется процессом скользящего среднего порядка m.

Заметим, что веса 1, 1, 2,..., m, на которые умножаются, не обязаны давать в сумме единицу или хотя бы быть положительными [2].

Если определить оператор скользящего среднего порядка m как ( B ) = 1 1B 2 B 2 K m B m, то модель скользящего среднего можно сжато записать, как xi = ( B ) i. Она содержит m + 2 неизвестных параметра:

µ, 1,..., m, 2, которые должны на практике оцениваться по наблюдениям.

Для достижения большей гибкости в подгонке моделей к наблюдаемым временным рядам иногда целесообразно объединить в одной модели и авторегрессию, и скользящее среднее. Это приводит к комбинированной модели авторегрессии - скользящего среднего [2] ( B ) xi = ( B ) i, в которой имеется p + q + 2 неизвестных параметра:

или µ; 1,K, p ; 1,K, q ; 2, оцениваемых по наблюдениям.

На практике часто оказывается, что адекватное описание наблюдаемых временных рядов достигается при помощи моделей авторегрессии, скользящего среднего или комбинированной модели, в которых p и q не больше, а часто и меньше 2 [4].

3.5. Методы моделирования случайных процессов Большой класс СП, имеющих место в информационно-измерительных системах, системах автоматического управления, а также в каналах связи, подверженных воздействию случайных возмущений описывается с помощью ДУ вида векторный стационарный СП.

При исследовании на ЭВМ системы (3.25) необходимо получать реализации СП. Методы моделирования СВ рассматривались в главе 2. Ниже приводятся некоторые распространенные на практике методы моделирования гауссовских стационарных СП: метод формирующего фильтра (п. 3.5.1), метод скользящего суммирования (п. 3.5.2) и рекуррентные моделирующие алгоритмы (п. 3.5.3).

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

Для линейных стационарных систем, находящихся под воздействием гауссовских стационарных случайных возмущений, могут быть получены алгоритмы моделирования, лишенные методических ошибок. Эти алгоритмы рассматриваются в п. 3.5.3. Они основаны на методе дискретизации линейных стохастических уравнений. Метод дискретизации дает сравнительно простые и легко реализуемые алгоритмы моделирования гауссовских векторных и скалярных СП с дробно-рациональным спектром высокого порядка. В п. 3.5. метод дискретизации применяется для процессов с типовыми КФ. Помимо задач цифрового моделирования алгоритмы дискретизации оказываются полезными при расчетах корреляционных характеристик линейных систем и применении методов оптимальной фильтрации к обработке СП.

В настоящее время разработан ряд методов моделирования гауссовских стационарных СП x ( t ) с заданными характеристиками: математическим ожиданием mx, КФ Rx ( ) или спектральной плотностью S x ( ). При решении задач моделирования в целях удобства зачастую считают математическое ожидание нулевым, а дисперсию 2 - единичной.

Использование преобразования позволяет получить процессы с требуемыми значениями этих характеристик.

S x ( ) описывает полностью рассматриваемый класс процессов. Обе характеристики связаны взаимно однозначно преобразованиями Фурье [36]:

Обычно задача формулируется следующим образом. По известным характеристикам процесса (математическому ожиданию, дисперсии и КФ или спектральной плотности) требуется построить вычислительный алгоритм, позволяющий получать на ЭВМ реализации СП x ( t ) или последовательностей xk, k = 0, 1, 2,.... В гауссовском случае модель процесса, заданная математическим ожиданием и КФ, является полностью определенной. Шаг дискретизации t может быть не равен шагу интегрирования h системы (3.25).

Известные методы можно разбить на две большие группы: точные (метод рекуррентных алгоритмов дискретизации) и приближенные (методы формирующего фильтра, скользящего суммирования). В точных методах отсутствует методическая ошибка по КФ, т. е. КФ Rx [l ] = M [ xk +l xk ] последовательности xk равна дискретным значениям Rx ( k t ) КФ Rx ( ) моделируемого процесса с непрерывным временем [41].

Для приближенных методов равенство заданных и воспроизводимых на ЭВМ характеристик выдерживается не точно, с некоторой погрешностью [41].

В настоящее время практически отсутствуют работы по анализу погрешностей приближенных методов моделирования, поэтому основным и наиболее надежным способом контроля приближенных алгоритмов является статистическая обработка моделируемых реализаций [7, 28, 35, 39].

3.5.1. Метод формирующего фильтра преобразующую СП ( t ) вида белого шума в СП x ( t ) с заданными статистическими характеристиками [1, 41]. Белый шум представляет собой стационарный СП с постоянной спектральной плотностью S0. Его КФ имеет вид где ( ) - дельта-функция Дирака, определяемая соотношениями Полагаем процесс ( t ) - гауссовским, нормированным условием S0 = 1 2, M ( t ) = 0. Чтобы найти передаточную функцию формирующего фильтра ( p ), спектральную плотность процесса x ( t ) представляют в виде произведения двух комплексно сопряженных сомножителей:

Формирующий фильтр с передаточной функцией ( p ) = S ( p ) S должен быть устойчивым. Отметим, что моделируемый процесс является стационарным с заданной спектральной плотностью лишь при t. Для дробно-рациональной спектральной плотности S x ( ) функция ( p ) имеет вид где Fm ( p ), H n ( p ) - полиномы степени m, n, m < n. Ей соответствует ДУ, записанное в операторной форме:

От этого уравнения с помощью известных преобразований легко перейти к системе ДУ первого порядка [41]. Стационарный СП x ( t ) может быть представлен первой компонентой n -мерного марковского процесса x ( t ) = x1 ( t ),..., xn ( t ), x1 ( t ) = x ( t ), удовлетворяющего уравнению матрица A и вектор B равны где bnm = 0, Процесс x ( t ) - гауссовский с нулевым средним. После окончания переходного процесса в уравнении (3.31) корреляционная матрица = M xxT, установившегося стационарного процесса находится из уравнения:

Белый шум с бесконечно большой дисперсией является абстрактным физически нереализуемым процессом. Для моделирования формирующего фильтра на ЭВМ разработаны различные способы, требующие предварительных вычислений [41].

Приведем приближенный и достаточно простой метод интегрирования на цифровой ЭВМ уравнений формирующего фильтра [7, 41]. На ЭВМ моделируется дискретный белый шум k ~ N ( 0, 1), k = 0, 1, 2,..., k некоррелированы. Рассмотрим ступенчатый процесс h ( t ) с шагом h, порождаемый дискретным белым шумом Спектральная плотность процесса h ( t ) равна При h 0 и фиксированном диапазоне частот [ 0, 0 ] функция Sh ( ) стремится к постоянной спектральной плотности, причем максимальное по отклонение достигается на конце промежутка при = 0. Относительная погрешность в имитации процессом h ( t ) свойств белого шума характеризуется величиной где * - заданное значение погрешности. Из неравенства (3.35) получаем неравенство h h* = 2 3* 0, позволяющее выбрать величину h.

Величина 0 определяет тот частотный диапазон, в пределах которого необходимо воспроизводить спектральную плотность S x ( ) моделируемого процесса. Значение 0 находится из условия где >0 - заданная малая величина. При этом диапазон частот [ 0, 0 ] должен перекрывать полосу пропускания системы, на вход которой подается процесс x (t ).

Уравнение формирующего фильтра при моделировании на ЭВМ, получается из формулы (3.30) при S0 = h 2 и имеет вид где h - шаг интегрирования ДУ формирующего фильтра. При моделировании на ЭВМ от этого уравнения следует перейти к системе ДУ первого порядка.

Векторному уравнению (3.31) соответствует уравнение Начальные условия задаются такими, чтобы можно было исключить переходный процесс: x ( 0 ) ~ N ( 0, ), - корреляционная матрица, определяемая из уравнения (3.34). При нулевых начальных условиях следует отбросить начальный отрезок реализации длиной T0. Для процессов с типовыми нормированными КФ x ( ) (см. примеры в п. 3.5.3 и 3.5.4) приведены расчеты, показывающие, что величина T0 приближенно равна T0 (1... 3) k, где k - интервал корреляции процесса. Величина k определяется условием x ( k ) = 0.05, где в случае неоднозначности в качестве k берется наибольший из корней уравнения [41].

Примеры применения метода формирующего фильтра для процессов с типовыми характеристиками приведены в п. 3.5.4. При малых h КФ процессов (3.31) и (3.37), а также (3.30) и (3.36) приближенно равны. При h методическая ошибка стремится к нулю. Отметим, что при интегрировании системы ДУ методом Рунге-Кутта на каждом шаге несколько раз вычисляются правые части. Составляя программы, следует предусмотреть, чтобы при вычислении правых частей использовалась одна и та же для данного шага интегрирования случайная величина k. Приведенный метод, основанный на моделировании посредством цифровой ЭВМ формирующего фильтра, обладает методической ошибкой, величина которой уменьшается при h 0.

3.5.2. Метод скользящего суммирования Дискретные значения моделируемого процесса формируются в виде скользящей суммы с весовыми коэффициентами a j. Существует ряд способов определения a j [1, 41]. Один из них основан на применении интеграла свертки где ( t ) - нормированный белый шум; g ( ) - весовая функция формирующего фильтра. Функция g ( t ) определяется формулой Формирующий фильтр с весовой функцией (3.40) имеет вещественную частотную характеристику функция (3.40) четна, поэтому непрерывный линейный фильтр с такой весовой функцией физически не реализуем. Однако это свойство не является препятствием для цифрового моделирования. Дискретизация интеграла (3.39) с шагом t дает следующие значения весовых коэффициентов:

Значения g ( j t ) вычисляются, как правило, с помощью численных методов [37]. При этом бесконечный верхний предел интегрирования в формуле (3.40) заменяют на конечный. Генерируемая последовательность xk имеет КФ Rx [l ] = M [ xk +l xk ], равную Истинная КФ имеет вид Функция Rx [l ] является интегральной суммой для интеграла (3.42). При условиях t 0, l t = = const КФ последовательности xk стремится к требуемой Rx ( ). Контроль правильности вычисления a j и выбора числа членов M (теоретически должно быть M = ) осуществляется путем расчета по формуле (3.41) функции Rx [l ] и сравнением ее с требуемой КФ. Поскольку последовательность (3.38) является гауссовской, то близость функций Rx [l ] и Rx ( ) означает близость заданного и моделируемого процессов на уровне конечномерных распределений. Метод скользящего суммирования пригоден для моделирования гауссовских процессов с произвольными спектральными плотностями.

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

Последовательное соединение исходной системы и формирующего фильтра образует эквивалентную систему [41].

Такая система может быть описана векторным ДУ где x - n -мерный вектор состояния, A = aij, B = bij - матрицы постоянных коэффициентов размеров n n и n m ; ( t ) - m -мерный гауссовский белый шум с нулевым средним и матричной КФ R ( ) = ( ) I m ; I m - единичная матрица. Начальные условия представляют собой гауссовский вектор x0 ~ N ( 0, 0 ), 0 - корреляционная матрица, x0 и ( t ) - независимые случайные векторы.

Дискретизация системы. Для моделирования системы (3.43) на цифровой ЭВМ перейдем к дискретной модели. Выразим x ( t ) через матричный экспоненциал Используя свойство e t + t в следующем виде:

где t - шаг дискретизации.

Соотношению (3.46) соответствует рекуррентное уравнение Векторы k +1 являются независимыми гауссовскими k +1 ~ N ( 0, R ) и не зависят от x k. Корреляционная матрица R = M k +1T +1 имеет вид Представим R в виде где B = bij - матрица порядка n r, r - ранг матрицы R.

Здесь матрица B находится по R с помощью алгоритма (2.20)-(2.22).

Тогда k +1 может быть выражен через k +1 ~ N ( 0, I r ) равенством k +1 = B k +1.

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

Алгоритм позволяет воспроизводить на ЭВМ случайные последовательности с заданными КФ [41].

Корреляционные моменты. Определение начальных условий для k = M {x k x k }, как следует из уравнения (3.51), определяется рекуррентным соотношением Примем далее, что матрица A в системе (3.43) – матрица Гурвица, т. е. ее собственные числа удовлетворяют условию Re j < 0, j = 1,2,..., n.

Отсюда, учитывая некоррелированность процесса ( s ) и его независимость от x0, после взятия математического ожидания получаем формулу Rx ( t, t ) = M {x ( t ) x ( t )} процесса x ( t ). При t = k t справедливо равенство После замены переменной интегрирования t s = u и функция примет вид Для матриц Гурвица имеем последовательно равенства:

Эти равенства означают, что в системе (3.43) асимптотически устанавливаются стационарные СП с нулевым средним и корреляционной Матрица находится предельным переходом в уравнении (3.52) при k. Переходя в формуле (3.43) к пределу, получаем линейное алгебраическое уравнение Матрицу можно также найти из соотношения (3.34). Если положить Это означает, что корреляционные свойства процесса x ( t ) не изменяются с течением времени. В системе (3.43) отсутствуют переходные процессы, а стационарный процесс устанавливается, начиная с момента t = 0, поэтому для устранения переходных процессов в уравнении (3.51) необходимо положить Матрица A = e A t находится из матричной системы ДУ интегрированием выражения (3.54) на промежутке [ 0, t ] при начальном условии X ( 0 ) = I n. После интегрирования получаем X ( t ) = A. Для приближенных вычислений можно пользоваться формулой (3.44), ограничиваясь конечным числом членов:

Рассмотрим, какой вид имеет уравнение (3.51) при малом шаге t. В (3.55) отбросим члены второго и высшего порядка малости. Тогда получим e At I n + A t, R BBT t. После подстановки этих выражений в формулу (3.51) следует приближенное равенство При стремлении t к нулю получаем Это выражение соответствует формуле (3.37), выведенной ранее другим способом.

Моделирующий алгоритм. Алгоритм цифрового моделирования включает следующие операции [41]:

1. Модель динамической системы или формирующего фильтра приводится к виду (3.43).

2. Вычисляется n -мерная матрица A. Для этого n раз интегрируется система где e j - j -й столбец единичной матрицы I n, x - n -мерный вектор. Решение системы (3.56) в момент t дает j -й столбец матрицы A.

3. Вычисляется матрица R = Rij. В соответствии с формулой (3.49) ее элементы равны где wil ( ) - компоненты матрицы w ( ).

Для определения Rij ) m раз интегрируются системы ДУ:

Здесь Bl - l -й столбец матрицы B, xi ( ) = wil ( ). В результате суммирования величин Rij ) ( t ) по l находятся элементы симметрической матрицы R.

4. По рекуррентному алгоритму (3.52) вычисляется корреляционная матрица. В качестве 0 можно принять любую неотрицательно определенную симметрическую матрицу, например нулевую. Выбор оказывает влияние только на время переходного процесса. Итеративный процесс (3.52) заканчивается, когда матрица k примет с заданной точностью установившееся значение.

5. С помощью алгоритма (2.20)-(2.22) вычисляется матрица B размера n r, определяемая соотношением (3.50).

6. Разыгрывается вектор x0 начальных условий; x0 ~ N ( 0, 0 ) или x0 ~ N ( 0, ) при моделировании стационарных процессов и установившихся режимов движения (3.43). Для этого также применяются алгоритмы (2.20)Разыгрыванием начального условия заканчивается подготовительный этап вычислений.

7. Выполняется цифровое моделирование динамической системы или стационарного СП по формуле (3.51).

Сравнивая алгоритм (3.51) с методами моделирования стационарных процессов, рассмотренными в п. 3.5.1 и п. 3.5.2, можно отметить следующее:

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

при моделировании отсутствуют переходные нестационарные процессы;

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

С использованием дискретной модели (3.51) могут решаться задачи статистического анализа линейных стационарных систем. Дисперсии и корреляционные моменты фазовых координат xi системы (3.43) являются в установившемся режиме. Эти матрицы находим из уравнений (3.52) и (3.53).

Рассмотрим матрицу Ее элементы представляют собой КФ и взаимные КФ. Как следует из уравнения (3.51) В установившемся режиме для стационарных процессов соотношение (3.59) принимает вид При большом n в выражении (3.43) применение формул (3.52), (3.53) и (3.59), (3.60) позволяет значительно сократить объем промежуточных вычислений по сравнению с известными методами анализа непрерывных систем [23, 32].

Пример 3.1. Частным случаем системы (3.43) является марковский гауссовский стационарный процесс первого порядка [5, 41]. Его нормированная КФ равна ( ) = e, уравнение (3.43) имеет вид Формулы (3.47)-(3.50) дают для коэффициентов дискретной модели (3.51) следующие выражения:

Модель (3.51) примет вид Пример 3.2. Рассмотрим стационарный СП x ( t ), имеющий спектр второго порядка [41] Используя соотношения (3.32) и (3.33), представим его в виде компоненты марковского процесса где b1 = x0, в уравнении (3.63) необходимо, чтобы начальное условие моделировалось в виде двумерного гауссовского вектора x0 ~ N ( 0, ), где корреляционная матрица = ij находится из (3.34) или (3.53). Решение уравнения (3.34) для данной системы дает следующие значения ij :

Для построения алгоритма (3.51) найдем матрицы A = e At, R и B. ДУ, необходимые для определения A ( ) = a ij ( ) имеют вид Первый столбец матрицы A ( ) получаем, решая эту систему при начальных условиях x1 ( 0 ) = 0, x2 ( 0 ) = 0, второй - при начальных условиях x1 ( 0 ) = 0, x2 ( 0 ) = 1. В результате имеем:

Из формул (3.49) и (3.57) следует, что матрица R равна где W ( ) - решение системы (3.65) с начальными условиями x1 ( 0 ) = b1, x2 ( 0 ) = b2. В силу линейности системы (3.65) имеем Отсюда получаем формулы для элементов R :

Матрицу B = b определяем из соотношения (3.50) с помощью формул (2.20)-(2.22). Элементы B равны:

В результате окончательно получаем Приведенные в данном примере формулы позволяют моделировать СП с различными типами КФ, встречающиеся на практике [4, 41].

3.5.4. Моделирование стационарных процессов с типовыми корреляционными функциями Во многих радиотехнических системах наблюдаются сигналы, которые достаточно хорошо описываются моделями стационарных СП с типовыми КФ [4]. Возмущения в динамических системах часто задаются также в виде гауссовских стационарных процессов с типовыми КФ и дробно-рациональными спектральными плотностями [1, 41].

Основываясь на рассмотренных в пп. 3.5.1-3.5.3 методах, получим алгоритмы для моделирования гауссовских стационарных СП с некоторыми типами КФ. В табл. 3.1 приведены КФ Ri ( ), спектральные плотности Si ( ) СП и соответствующие им передаточные функции i ( p ) формирующих фильтров.

Приведенные в табл. 3.1 типовые КФ имеют следующие случайные возмущения, встречающиеся в приложениях: атмосферная турбулентность;

шумы/помехи в следящих системах и информационно-измерительных устройствах; неоднородности земной поверхности; сейсмические нагрузки;

характеристики потоков событий и др. [41] Системы ДУ вида (3.37) для моделирования формирующих фильтров на ЭВМ представлены в табл. 3.2. При этом x1 ( t ) = x ( t ). Там же даны значения шага интегрирования h*, при котором обеспечивается заданная точность воспроизведения спектральной плотности (3.35) в установленном диапазоне частот. КФ R5 ( ) = 2 e (1 + ) получается в результате предельного перехода при 0 из КФ R3 ( ).

Моделирующие алгоритмы получаются из формул табл. 3.2 при = 0. Для исключения переходного процесса начальные условия в ДУ (табл. 3.2) следует задавать как реализацию случайного вектора x0 ~ N ( 0, ). Подставляя в (3.64) параметры КФ, определяем корреляционные моменты ij ; 11 = 2, значения 12 и 22 приведены в табл. 3.3.

Получим дискретные модели, позволяющие моделировать процессы с типовыми КФ без методических ошибок. Процесс с экспоненциальной КФ приведенными в табл. 3.1, имеют спектры второго порядка вида (3.62).

Представление в виде компоненты марковского процесса имеет вид (3.63), где 1 = 2, а значения остальных параметров приведены в табл. 2.4. Дискретная модель определяется уравнениями (3.69), где значения aij, bij определяются подстановкой данных табл. 3.4 в формулы (3.66)-(3.67). Окончательные выражения для aij, Rij через параметры КФ приведены в табл. 3.5.

Коэффициенты bij связаны с Rij формулами (3.68). Стационарные СП с дробно-рациональной спектральной плотностью (3.62) могут моделироваться с помощью уравнения типа авторегрессии - скользящего среднего [2, 41]:

Для определения параметров уравнения (3.70) найдем спектральную плотность последовательности k :

плотности S ( ) = Sij ( ) последовательности (3.69) и определяется формулой Обратная к матрице G ( z ) порядка 2 2 легко находится с помощью присоединенной матрицы. Подставив в формулу (3.71) элементы матриц A и R, получим выражение для спектральной плотности:

где в знаменателе – определитель матрицы G ( z ), постоянные ci, ai равны Для процессов с типовыми КФ Ri ( ), i = 2, 3, 4, отсюда получаем Здесь для КФ R4 ( ) нужно положить = 0. Для определения коэффициентов b1, b2 выполним факторизацию числителя, т. е. представим его в виде Корни трехчлена c0 z 2 + c1 z + c0 имеют вид v1, 2 = 1 ± ( c1 2c0 ) 1.

В качестве v1 примем тот из корней, который по модулю меньше единицы, он определяется формулой Второй корень равен v2 = 1 v1. Разлагая трехчлен на линейные множители, преобразуем P ( z ) к виду где c = c0 v1 - положительная величина. Отсюда следует Для того, чтобы избавиться от переходного процесса, необходимо разыгрывать начальные условия ( x0, x1, 2, 1 ) как гауссовский четырехмерный случайный вектор с нулевым средним и корреляционной матрицей Для удобства пользования дискретный алгоритм моделирования процессов с типовыми КФ Ri ( ) и его параметры представлены в табл. 3.6. Приведенные в таблице выражения совпадают с формулами табл. 3.2.

3.6. ИДЕНТИФИКАЦИЯ И АНАЛИЗ АДЕКВАТНОСТИ

АВТОРЕГРЕССИОННЫХ МОДЕЛЕЙ СЛУЧАЙНЫХ

ПРОЦЕССОВ

Рассмотрим построение статистической модели СП путем определения параметров модели. Типичными моделями СП, которые могут легко реализованы на ЭВМ являются модели авторегрессии-скользящего среднего, рассмотренные в п. 3.4. Данная задача также называется идентификацией модели и решается она путем статистического оценивания параметров СП [2, 15].

В частности, идентификация модели осуществляется при исследовании статистических характеристик радиопомех и замираний в каналах связи, представимых в виде последовательности дискретных отсчетов x1,..., xn, поскольку их огибающие удобно представлять в виде временных рядов с заданными корреляционными свойствами [3, 4].

На основе полученной выборки отсчетов x1,..., xn выполняется идентификация модели СП путем оценивания параметров модели, а именно дисперсии и коэффициентов корреляции:

После этого составляется система уравнений Юла-Уокера (см. п. 3.4.1), которая, например, для АР модели 2-го порядка будут иметь вид В систему уравнений (3.73) подставляются коэффициенты корреляции из (3.72) и находятся соответствующие коэффициенты авторегрессии. При этом в случае авторегрессии 1-го порядка вместо системы будет лишь равенство 1 = R (1).

Диагностическая проверка адекватности модели, заключается в выражении белого шума через коэффициенты авторегрессии, полученные из (3.73) и входные наблюдения zt с помощью исходных уравнений авторегрессии. При диагностике АР моделей в качестве критерия адекватности принимается критерий проверки нулевой гипотезы относительно независимости соседних отсчетов разностей i = xi xi 1 2 xi 2. Для этого используется следующий тест где n - количество отсчетов в выборке, r - выборочный коэффициент корреляции остатков i модели авторегрессии. При этом имеем Обычно анализируется выборка размером n > 200, уровень значимости принимается равным 0.05. По таблице критических точек распределения Стьюдента по заданному уровню значимости и числу степеней свободы k = n 2 выбирается критическая точка tкр =1.96 для двухсторонней критической области. Если о независимости остатков, т. е. об адекватности модели. В противном случае требуется вычислить еще один коэффициент корреляции исходной выборки R ( k ) и решить систему Юла-Уокера более высокого порядка, после чего повторить всю процедуру проверки адекватности. Очевидно, что чем выше порядок АР модели, тем она более гибкая в смысле аппроксимации параметров исходной выборки данных.

3.7. Заключительные замечания Таким образом, задача цифрового моделирования СП формулируется как задача нахождения алгоритмов (по возможности наиболее простых), позволяющих получать на ЭВМ дискретные реализации (выборочные функции) моделируемых процессов. Это самостоятельная и довольно сложная задача синтеза дискретных СП, имитирующих непрерывные процессы с заданными статистическими характеристиками. Она решается путем отыскания удобных для реализации на ЭВМ линейных и нелинейных преобразований, с помощью которых можно превратить независимые равномерно или нормально распределенные случайные числа, вырабатываемые датчиком случайных чисел, в случайные последовательности с требуемыми вероятностными свойствами [1, 4, 7, 28].

Ввиду ограниченности объема данного пособия за рамками изложения остались многие специальные методы моделирования СП, такие как метод канонических разложений [1, 6], одномерные и двумерные модели разрывных СП [16], методы имитации импульсных помех на ЭВМ [4, 16] и др.

4. МОДЕЛИРОВАНИЕ СЛУЧАЙНЫХ ПОЛЕЙ

Решение различных задач математического моделирования, возникающих при проектировании телекоммуникационных систем, а также во многих других областях, приводит к необходимости получения на ЭВМ реализаций многомерных массивов СВ с заданными вероятностными характеристиками – случайных полей. Модели случайных полей используются для имитации фона в оптико-электронных системах, при описании свойств рельефа земной поверхности, для моделирования метеополей, применяются в устройствах фильтрации и кодирования изображений. Типичными примерами применения моделей многомерных случайных полей в системах связи являются: описание фазовых искажений на двумерной сетке время-частота; описание пространственно-временных сигналов в системах связи с несколькими антеннами (MIMO – Multi-Input Multi-Output).

При решении задач обработки случайных полей важным этапом является выбор адекватной модели наблюдений [5, 12, 31]. В настоящее время не существует универсального способа формирования случайных полей с произвольно заданными характеристиками. Поэтому известные модели случайных полей соответствуют реальным случайным полям лишь по ограниченному числу параметров (форма КФ, распределение амплитуд и т.п.) [8, 11, 12, 31, 41]. Ниже рассмотрен ряд известных моделей, которые могут быть использованы для приближенного описания случайных полей при синтезе различных алгоритмов обработки.

Наиболее изученными являются АР модели случайных полей [5, 20, 31].

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

В настоящем пособии при рассмотрении моделей случайных полей используется идея представления спектрально-корреляционных характеристик многомерных сигналов в разделимой по пространственным координатам форме [31]. Это дает возможность использования достижений теории одномерных сигналов, приводя к относительно простым ММ, уже позволившим получить решение ряда задач статистической обработки данных [20].

Эффективным методом решения разнообразных задач обработки сигналов служит спектральный анализ [8, 31, 43]. К сожалению, существует лишь узкий класс так называемых «разделимых» случайных полей на многомерных сетках [31], для которых можно получить полезные для приложений аналитические соотношения. В частности, важнейший класс изотропных случайных полей дискретного аргумента не удается исследовать известными методами спектрального анализа. Это объясняется несоответствием декартовой системы координат в пространстве R0N точек с целочисленными координатами и естественной для изотропных случайных полей сферической системой координат в R N.

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

Задачи моделирования случайных полей относятся к сравнительно новой и еще мало изученной области в теории и практике статистического моделирования. Особую важность здесь приобретает проблема сокращения вычислительных затрат. Объемы вычислений и машинное время ЭВМ резко возрастают с ростом размерности случайного поля, а также с уменьшением шага дискретизации. Модели многомерных случайных полей (скалярных и векторных, гауссовских и негауссовских) теоретически мало изучены [41].

4.1. Моделирование многомерных дискретных случайных информационные комплексы, включающие в себя пространственные системы датчиков и цифровую вычислительную технику. Поэтому будем в основном рассматривать случайные поля с дискретными пространственными и временными переменными. Не ограничивая общности, будем считать, что случайные поля заданы на многомерных прямоугольных сетках с единичным шагом. На рис. 4.1, а и 4.1, б изображены двумерная и трехмерная сетки. В общем случае случайное поле задано в узлах n -мерной сетки В зависимости от физической природы значения отсчетов случайного поля могут быть скалярными (например, яркость монохроматического изображения) или векторными (поле скоростей, цветные изображения, поле смещений). Если обозначить через x j значение отсчета случайного поля в узле (пикселе) j, то случайное поле есть совокупность этих значений на сетке: X = {x j : j }.

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

Рис. 4.1. Двумерная (а) и трёхмерная (б) прямоугольные сетки Если xi j - некоторая физическая величина в точке ( i, j ), то совокупность являющегося случайной функцией x ( i, j, ) двух пространственных переменных и СВ. Рассматривая подобные процессы развивающимися во времени, к пространственным переменным ( i, j ) нужно добавить еще и временную переменную i. При этом поле X = { xii j } становится трехмерным.

четырехмерного поля и т.д.

Если требуется временную переменную выделить особо, то будем ее записывать сверху: X = { x j : j, i I }. Такое случайное поле задано на прямом произведении I сеток и I, где I – множество значений временного индекса. Сечение x = {x j : j }, т. е. совокупность отсчетов случайного поля при фиксированном значении временного индекса i, называется i -м кадром случайного поля. Каждый кадр задан на сетке.

Например, на рис. 4.1, б изображено три двумерных кадра.

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

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

Во-вторых, значительно усложняется понятие марковости.

СП можно представить развивающимся во времени, математическим выражением такого развития и является модель (3.1). Для марковских последовательностей временной интервал может быть разбит любой точкой i на + = { x k : k > i}. Однако случайное поле определено на n-мерной области, для геометрического разбиения которого на две части и требуется, по меньшей мере, ( n 1) -мерная область. Свойство марковости случайного поля состоит в том, что для любого множества (из некоторого класса множеств) СВ, входящие в, условно независимы от СВ, входящих в, при известных значениях. Назвать, и прошлым, настоящим и будущим можно весьма условно. Тем не менее марковское свойство позволяет представить случайное поле также формирующимся во времени от через к, причём с течением времени перемещается по. Например, если в качестве брать строки двумерной сетки, то поле X можно представить формирующимся построчно [20].

Дальнейшее развитие этой идеи позволяет обобщить АР модели случайных последовательностей на случайные поля. В связи с этим рассмотрим в общем виде задачу рекуррентного формирования случайного поля { x j ; j = ( j1 j2... jN ) J } на N -мерной прямоугольной сетке J = { j, jk = 1 + M k, k = 1, 2,..., N }. При этом предполагается, во-первых, существование некоторой процедуры последовательного перебора точек j J, т.е. правила линейного упорядочения точек j, l J, на основе которого можно сказать, что элемент j предшествует элементу l ( j < l ) или наоборот. Во-вторых, должен быть задан алгоритм, определяющий, каким образом очередное значение случайного поля x j может быть найдено на область индексов l J, предшествующих очередному элементу j. Такую область G j конечных размеров обычно называют областью локальных состояний [8, 12, 31]. Наконец, для формирования случайного поля j с определенными вероятностными характеристиками на каждом шаге рекуррентных вычислений функция x j = j ( xl, l G j ) должна включать в Таким образом, представление случайного поля на основе рекуррентной процедуры должно иметь следующий вид где G j - области элементов l J, на которых уже определены предыдущие значения случайного поля { xl } ; j ( xl ; l ), j J, вообще говоря, нелинейные скалярные или векторные функции двух тензорных аргументов. Наиболее простым частным случаем (4.1) является линейное стохастическое уравнение с белым гауссовским СП { l }, соответствующее известному уравнению авторегрессии - скользящего среднего [2] для СП. Однако в отличие от своего одномерного аналога, свойства случайного поля x j, j J, порождаемого (4.2), в настоящее время изучены не полностью даже для моделей (4.2) с постоянными коэффициентами j, l = l, j, l = l и не изменяющимся видом областей G j = G и Y j = Y :

Если порядок формирования случайной последовательности x 0, x1, x 2,...

обычно соответствует наблюдаемым во времени значениям, то порядок формирования случайного поля X = { x j : j } требует дополнительного определения. Для этого нужно линейно упорядочить узлы сетки, тогда про любые два элемента поля можно сказать, что один из них предшествует другому. Если xi предшествует x j, то будем отмечать это как ( i ) < ( j ), т. е. номер элемента xi меньше номера x j при данной развертке. Существует множество вариантов такого упорядочения. В двумерном случае чаще всего применяются пилообразная и треугольная развертки (рис. 4.2).

Рис. 4.2. Пилообразная (а) и треугольная (б) развертки В результате развертки поле преобразуется в случайную последовательность. Предположим, что она является марковской порядка s, т. е.

условная ПРВ любого xi относительно всех предшествующих ему элементов зависит только от элементов некоторого конечного отрезка i = {x j : (i ) s ( j) < (i )}. Множество i называется глобальным состоянием.

В двумерном случае оно при пилообразной (и треугольной) развертке включает в себя несколько последних строк и показано на рис. 4.3. Следовательно, можно представить xi в каузальном виде как функцию элементов глобального состояния и возмущения i :

Полученное выражение представляет АР модель случайного поля. Однако использовать (4.4) для представления полей на сетках больших размеров трудно, а для бесконечных сеток – невозможно ввиду большого или даже бесконечного числа аргументов функций i.

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

Рис. 4.3. Области локального и глобального состояний В результате поле X может быть представлено АР моделью которая во многих случаях может быть приемлема для решения прикладных задач. Конечно, может оказаться, что даже область локального состояния Li слишком велика, и возникают значительные технические трудности при имитации или обработке полей. В таких ситуациях можно Li уменьшить до приемлемых размеров, используя полученную модель (4.5) как некоторое приближение к реальным физическим объектам.

4.1.2. Авторегрессионные модели случайных полей Рассмотрим линейную гауссовскую АР модель коэффициенты;

стандартных гауссовских СВ. Одним из первых подобную модель применительно к оценке плоских изображений исследовал Хабиби [20]:

При этом первый элемент поля формируется как x1,1 = X 1,1 ;

первый столбец xi1,1 = 1 xi1 1,1 + X 1 1 i1,1 ;

первая строка x1, i2 = 2 x1, i2 1 + X 1 2 1, i2.

Схема вычислений для этой модели представлена на рис. 4.4, а.

Порождаемое поле имеет разделимую экспоненциальную КФ Рис. 4.4. Схемы формирования авторегрессионных СП Перепишем (4.7) в виде где z11 xi1,i2 = xi1 1,i2, z2 1 xi1,i2 = xi1,i2 1 - операторы сдвига вдоль соответствующих направлений.

Для n -мерного случая будем иметь или Порождаемое поле является анизотропным и имеет множительную КФ На основе модели (4.7) разработано большое количество алгоритмов фильтрации случайных полей. Однако она, как и ее многомерный вариант (4.9), имеет существенный недостаток – множительность (факторизуемость) КФ. В двумерном случае элементы поля, одинаково коррелированные с элементом xi1i2, расположены на ромбе с центром в ( i1, i2 ), а в многомерном случае – на ромбоиде, хотя более естественными сечениями КФ для реальных полей были бы эллипс и эллипсоид (рис. 4.5). Для частичного скругления сечений КФ можно отказаться от частного вида весовых коэффициентов (4.9), а также расширить область локальных состояний. Однако такое расширение приводит к резкому увеличению числа слагаемых в (4.6).

Существенного скругления сечений КФ случайного поля позволяют достичь модели с кратными корнями характеристических уравнений одномерных АР моделей [45]. Решим сначала эту задачу для одномерной модели, для чего рассмотрим одномерную авторегрессию длины M :

с нулевым математическим ожиданием и дисперсией. Решение задачи синтеза будет заключаться в том, чтобы по заданным 1 – корню характеристического уравнения, n – его кратности и x – требуемой дисперсии поля определить неизвестные коэффициенты авторегрессии Модели (4.6) и её частному случаю (4.10) соответствует пространственный линейный фильтр с передаточной функцией где z j = z1 1 z 2 2 K z N N. При этом спектральная плотность случайного поля X записывается следующим образом:

КФ поля X может быть найдена с помощью обратного z -преобразования спектральной плотности:

где C N = { z = 1} – единичная полиокружность в многомерном комплексном пространстве.

В случае с кратными корнями уравнение (4.10) можно записать в операторной форме следующим образом:

где z 1 – оператор сдвига. Учитывая, что действие оператора сдвига на i -й элемент последовательности определяется как z k xi = xk i, перепишем (4.13) в явном виде:

Сравнение (4.13) и (4.14) дает возможность записать выражение для коэффициентов j = j (, n ) :

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

Одной из задач статистического анализа модели, является получение ее КФ [45]. Найдем вначале нормированную КФ, т. е. будем полагать 2 = = 1.

Для решения этой задачи воспользуемся формулой (4.12), записанной для одномерного случая Так как подынтегральная функция имеет в точке z = полюс порядка n, то интеграл представляет собой коэффициент c n ее разложения в ряд Лорана, и поэтому может быть найден с использованием методов теории вычетов Используя правила дифференцирования произведения функций, после предельного перехода получаем где а коэффициент = (,n ) находится из условия R(0 ) = 1:

Формулы (4.16)-(4.18) дают, при заданных и n, общий вид нормированной КФ одномерной модели (4.10). Для того чтобы получить КФ при не равных единице дисперсиях 2 и, необходимо домножить правую часть (4.18) на 2 /. Тем самым, получаем выражение для коэффициента :

Формулы (4.15) и (4.19) полностью определяют неизвестные коэффициенты одномерной АР модели (4.10) с кратными корнями характеристических уравнений [45].

Рассмотрим теперь случай N измерений. Модель случайного поля, при заданной дисперсии 2, полностью определяется вектором параметров ( 1, 2, K, N ) и вектором кратностей (n1, n2, K, n N ) характеристических корней.

Пусть многомерное разделимое случайное поле порождается следующими АР уравнениями, записанными в операторной форме:

где N – размерность поля; k и nk – параметр и кратность корней модели вдоль k -й оси; – сетка, на которой определено поле X.

Определим коэффициенты авторегрессии для многомерной модели с кратными корнями. Для этого раскроем в (4.20) скобки:

на N -мерном параллелепипеде размера (n1 + 1) (n 2 + 1) K (n N + 1). Из (4.20) и (4.21) следует, что коэффициенты j являются произведениями соответствующих коэффициентов kjk одномерных авторегрессий вдоль k -й оси:

где j = ( j1, j 2, K, j N ), j k = 1K n k. Коэффициент многомерной модели (4.21) находится аналогично:

где k – соответствующий нормированный одномерный коэффициент.

Таким образом, выражение (4.21) дает общий вид АР модели многомерного разделимого гауссовского случайного поля, а формулы (4.22) и (4.23) полностью определяют ее коэффициенты, т. е. задача синтеза модели решена. КФ модели (4.21), как уже отмечалось, является произведением КФ соответствующих одномерных авторегрессий:

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

Рассмотрим некоторые примеры. Хорошо изученная трехточечная модель Хабиби (4.7) является частным случаем 2-мерной модели (4.21) кратности (1,1), причем значение параметра задает коэффициент корреляции соседних элементов. Легко видеть, что ее КФ, вычисленная по формулам (4.16), (4.24), совпадает с (4.8):

Для 2-мерной модели кратности (2,2) КФ записывается более сложно [45]:

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

Это подтверждается также видом сечений равного уровня КФ, приведенных на рис. 4.5.

параметры (0.9, 0.9) параметры (0.9, 0.9) параметры (0.9, 0.9) Полезный во многих приложениях вариант модели многомерного случайного поля можно получить, взяв за основу одномерную авторегрессию с кратными корнями. Пусть необходимо сформировать реализации N -мерного = 1, M 1 1, M 2... 1, M N. Это можно сделать следующим образом:

x = x1, x 2,..., x M k ; k = 1...N на основе модели с кратными корнями (4.14) ( n = (n1, n2,..., n N ) – вектор кратностей вдоль соответствующих осей). Далее, элемент x j поля X получается перемножением соответствующих элементов одномерных последовательностей:

Здесь R k () – КФ (4.16) одномерной АР модели вдоль k -й оси. Таким образом, данная модель является разделимой моделью случайного поля, причем ее КФ совпадает с КФ N -мерной АР модели с кратными корнями. Закон распределения вероятностей такого случайного поля оказывается негауссовским и, в общем случае N, достаточно сложным. Вместе с тем, рассмотренная модель отличается простотой и очень малыми вычислительными затратами. Действительно, даже для простейшего случая N = 2 одномерная ПРВ случайного поля находится как ПРВ произведения двух гауссовских СВ в виде [38]:

где K 0 () – модифицированная функция Бесселя третьего рода нулевого порядка.

Проведем оценку вычислительной сложности предложенной модели. Для получения реализации N -мерного СП, определенного на сетке размером M 1 M 2 K M N элементов, требуется O ( nk + 1) M k операций умножения.

4.2. Алгоритмы моделирования непрерывных случайных полей Для специальных классов распределений задача моделирования случайного вектора может быть существенно упрощена. В приложениях часто приходится моделировать изотропные случайные векторы [4, 41]. К этому классу распределений приводит, например, моделирование изотропных случайных полей. Случайный вектор x = ( x1,..., xm ) называется изотропным, если его орт e = x x распределен равномерно на поверхности m -мерной сферы и не зависит от распределения величины x. Здесь x = x12 +... + xm - модуль вектора. Для того, чтобы случайный вектор x был изотропным, необходимо и достаточно, чтобы его ПРВ зависела только от модуля x, т. е.

где w ( x ) - некоторая неотрицательная функция, x ( 0, + ).

Переходя к сферической системе координат r, x, можно представить ПРВ изотропного случайного вектора в виде где S1 = 2 ( m 2 ) - площадь поверхности единичной сферы в m -мерном пространстве; r и x соответственно длина вектора и его проекция на сферу единичного радиуса; w0 ( r ) - плотность распределения x :

Моделирование изотропного случайного вектора осуществляется по формуле x = x e, где x - СВ с плотностью w0 ( r ), e = ( e1,..., em ) реализация изотропного направления в m -мерном пространстве R m. Вектор e является ортом вектора x, e = 1. При m = 2 изотропное направление задается на плоскости полярным углом, равномерно распределенным на [ 0, 2].

Моделирующий алгоритм имеет вид Рассмотрим двухпараметрическую модель случайного поля [41]. Пусть требуется получить в области D R m реализацию однородного вещественного случайного поля X ( x ) с заданным математическим ожиданием mX, дисперсией 2 и КФ K X ( x ) = 2 X ( x ). Здесь X ( x ) - нормированная КФ случайного поля. Следует отметить, что задача моделирования случайного поля по заданным двум первым моментам математически однозначно определена лишь для гауссовского поля. Для негауссовских полей данная задача не имеет однозначного решения [11].

где v R m - m -мерный случайный вектор с ПРВ равной S X ( u ) ; z - СВ с ПРВ детерминированный вектор; x - аргумент поля; vT ( x0 + x ) - скалярное произведение векторов v и x0 + x.

Моделирование случайного поля на основе параметрической модели (4.28) сводит задачу к получению реализаций m -мерного случайного вектора v, имеющего ПРВ равную спектральной плотности S X ( u ) моделируемого поля.

Моделирование вектора v значительно упрощается при имитации специальных классов полей. Рассмотрим распространенный на практике случай изотропных случайных полей, имеющих корреляционные и спектральные характеристики, зависящие от модулей соответствующих векторов:

При моделировании изотропных полей распределение вектора v также изотропно. Моделирование изотропного случайного вектора v сводится к моделированию его модуля v - СВ с плотностью распределения w0 ( u ) (4.26) и орта e, задающего изотропное направление в пространстве частот. Для двухи трехмерных полей плотности w0 ( u ) равны соответственно Алгоритмы моделирования v для двух- и трехмерных полей с типовыми характеристиками представлены в табл. 4.1. ПРВ v вычислялись по формулам (4.29) и указанным в табл. 4.1 типовым спектрам. Вывод моделирующих алгоритмов дан в главе 2. Для двумерных полей изотропное направление задается формулой (4.27); в трехмерном случае удобен алгоритм:

Примечание., 1, 2 ~ Uni [ 0, 1] ; 1, 2 ~ N ( 0, 1).

Моделирующие алгоритмы в многомерном случае можно найти в работе [13].

Как уже отмечалось, ряд моделей анизотропных полей описываются с помощью КФ множительного вида x = ( x1,..., xm ). В этом случае спектральная плотность равна и, следовательно, компоненты вектора v - независимые величины с ПРВ, последовательным разыгрыванием реализации его компонент v1,..., vm как независимых СВ.

В табл. 4.2 представлены алгоритмы моделирования компонент вектора v для полей, имеющих КФ вида (4.30) с типовыми характеристиками. Алгоритм моделирует распределение Коши, алгоритм 2 - нормальное распределение, алгоритм 3 - равномерное распределение на отрезке [, ]. Моделирующие алгоритмы 4 и 5 выведены в главе 2.

Номер алгоритма Примечание. ~ Uni [ 0, 1] ; ~ N ( 0, 1) ; = ±1 с вероятностью p = 0.5.

Ряд анизотропных случайных полей, в частности, поля, имеющие эллиптическую анизотропию, могут моделироваться как изотропные случайные поля. Анизотропия воспроизводится на ЭВМ изменением масштаба аргумента поля. Например, двумерное случайное поле с КФ вида заменой y ' = Ly сводится к изотропному случайному полю вида 1 (табл. 4.1), = 1 rx. Здесь ( x, y ) R 2, rx, ry - интервалы корреляции поля вдоль осей X, Y соответственно, L = rx ry - степень анизотропии. Обратной заменой y = y ' L, которая соответствует сжатию реализации вдоль оси Y (при L > 1 ) и растяжению (при L < 1 ), осуществляется переход от изотропного поля к модели (4.31).

Для поля со спектральной плотностью и КФ вида где B - положительно определенная матрица порядка m, реализации вектора v получаются моделированием m -мерного гауссовского вектора с нулевым математическим ожиданием и корреляционной матрицей B 1.

4.3. Заключительные замечания Рассмотренный круг моделей случайных полей не является исчерпывающим. Модели высших порядков могут обеспечить лучшую аппроксимацию реальных изображений. Однако выбор модели обусловлен, помимо требования точности описания, возможностью построения эффективного алгоритма обработки [8, 11, 31]. С учетом этого рассмотренные модели хорошо соответствуют реальным данным, обладают существенными преимуществами с вычислительной точки зрения, позволяя строить рекуррентные процедуры оценивания и использовать при обработке быстрые спектральные преобразования [8].

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

5. МОДЕЛИРОВАНИЕ СЛУЧАЙНЫХ ПОТОКОВ

И СИСТЕМ МАССОВОГО ОБСЛУЖИВАНИЯ

Системы массового обслуживания (СМО) отражают процессы функционирования многих современных сложных систем. Теория массового обслуживания получила свое первоначальное развитие благодаря запросам телефонии при изучении процесса загрузки телефонной станции [6, 18, 21].

Примеры СМО: телефонные станции, билетные кассы, системы организации транспорта (грузоперевозки, поток автомобилей через мост) и т. п.

Каждая такая система состоит из некоторого количества обслуживающих единиц: «каналов» обслуживания (например, линии связи). Работа любой СМО состоит в обработке (обслуживании) поступающего на нее потока требований или заявок. Заявки поступают одна за другой в случайные моменты времени.

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

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

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

5.1. Модели случайных потоков 5.1.1. Виды потоков и способы их задания В теории массового обслуживания вводится в рассмотрение модель потока событий. Потоком событий называется последовательность событий, следующих одно за другим в случайные моменты времени t1, t2,..., tk,...

(рис. 5.1). Потоки событий, происходящих в случайные моменты времени t 1, t 2 > t 1,..., t n > t n 1,..., являются специфичным классом СП. Они используются в качестве ММ при исследовании СМО, в задачах приема импульсных сигналов, в задачах надежности и др. [21] На рис. 5.1 изображена одна реализация потока событий. При этом в качестве событий могут быть, например, заявки на входе СМО и заявки на выходе СМО.

Потоки заявок могут быть однородными и неоднородными.

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

Возможны различные эквивалентные способы задания случайных потоков.

Наиболее общим способом представления характеристик потоков является задание многомерной ПРВ интервалов между моментами наступления событий w ( 1,..., n ), где k = tk tk1, t0 = 0.

При таком задании случайных потоков моделирование сводится к формированию на ЭВМ реализаций случайных чисел 1, 2,..., n с законом распределения w ( 1,..., n ).

Случайные потоки общего вида редко встречаются в приложениях.

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

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

В реальных СМО поток заявок обычно является случайным. Такие потоки характеризуются следующими тремя свойствами.

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

2. Поток событий называется потоком без последействия, если для любых, не перекрывающихся интервалов времени ( t1,t1 + t1 ) и ( t2,t2 + t2 ), где t1 > 0, t2 > 0, t2 t1 + t1, вероятность появления числа событий, попадающих на один из них, не зависит от числа событий, попадающих на другой интервал (рис. 5.3).

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

5.1.2. Простейший поток Если поток событий обладает тремя названными свойствами, то он называется простейшим (или стационарным пуассоновским) потоком.

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

Пусть имеется ряд независимых потоков 1, 2,..., n, имеющих одинаковое распределение. «Суммирование» потоков состоит в том, что все моменты появления событий переносятся на одну и ту же ось времени Ot.

Рассмотрим на оси Ot два неперекрывающихся отрезка (рис. 5.3). Каждая из точек, попадающих в эти отрезки, случайным образом может оказаться принадлежащей тому или иному потоку, и по мере увеличения n удельный вес точек, принадлежащих одному и тому же потоку (и, значит, зависимых), уменьшается, а остальные точки принадлежат разным потокам и появляются на отрезках независимо друг от друга. При увеличении n суммарный поток будет терять последействие и приближаться к простейшему. На практике достаточно к простейшему.

Важной характеристикой является интенсивность потока, определяемая как предел где 1 ( t0, t ) — вероятность того, что на интервале ( t0, t0 + t ) появятся одна или более заявок. Для стационарного потока его интенсивность не зависит от времени ( t ) = и равна среднему числу событий в единицу времени.

Рассмотрим на оси Ot простейший поток событий. Выделим произвольный участок времени длиной. Число точек m, попадающих на участок, распределено по закону Пуассона В частности, вероятность того, что не произойдет ни одного события P0 ( ) = exp ( ) ; вероятность появления ровно одного события P ( ) = e.

Важной характеристикой потока является закон распределения случайного интервала T между соседними событиями. Функция распределения F ( ) = P (T < ) для простейшего потока находится следующим образом:

Дифференцируя простейшего потока:

показательному закону m = M [ ] = 1, а дисперсия D = D [ ] = 1.

Следует заметить, что кроме характеристик входного потока заявок, режим работы СМО зависит еще от характеристик производительности самой системы: числа каналов n и быстродействия каждого канала. Одной из важнейших величин, связанных с системой, является время обслуживания одной заявки Tоб.

Рассмотрим СВ Tоб и обозначим Fоб ( t ) ее функцию распределения:

Fоб ( t ) = P (Tоб < t ), а wоб ( t ) = Fоб ( t ) - ПРВ. Для практики особый интерес представляет случай, когда величина Tоб имеет показательное распределение где параметр µ = 1 mtоб называется интенсивностью обслуживания.

Как известно параметр имеет смысл «плотности потока заявок».

Аналогично, величину µ можно характеризовать как «плотность потока освобождений» занятого канала. Представим себе канал, непрерывно занятый (бесперебойно снабжаемый заявками); тогда очевидно, в этом канале будет иметь место простейший поток освобождений с плотностью µ.

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

Потоки Пальма часто получаются в виде выходных потоков СМО. Если на какую-либо СМО поступает поток заявок, то он этой системой разделяется на два: поток обслуженных заявок и поток необслуженных заявок, который, в свою очередь, поступает на какую-либо другую СМО.

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

Потоки, у которых w1 ( ) = w ( ), определяются единственным законом распределения w ( ) и называются рекуррентными (стационарными) w2 ( ) = w3 ( ) =...wn ( ) = w ( ), называются рекуррентными потоками с запаздыванием. Они задаются двумя законами распределения w1 ( ) и w ( ).

Здесь wi ( ) характеризует ПРВ временного интервала между i -й и ( i + 1) -й заявками.

Рассмотрим моделирование потоков с ограниченным последействием. Для получения реализаций последовательности моментов наступления событий tk, k = 1, 2,..., достаточно сформировать последовательность реализаций k, k = 1, 2,..., СВ с заданными законами распределения wk ( ) соответственно и вычислить моменты наступления событий: tk = tk 1 + k.. Моделирование рекуррентных потоков упрощается еще и тем, что СВ k имеют одинаковый закон распределения.

Рис. 5.4. Получение потока Эрланга путем просеивания Пример потоков с ограниченным последействием – потоки Эрланга, которые образуются «просеиванием» простейшего потока (рис. 5.4). Если в простейшем потоке выбросить каждую вторую точку, то оставшиеся точки образуют поток, называемый потоком Эрланга первого порядка ( Э1 ). Такой поток является потоком Пальма, поскольку величины 1, 2,..., получаются суммированием независимых интервалов. Вообще, потоком Эрланга k -го порядка ( Эk ) называется поток, получаемый из простейшего, если сохранить каждую ( k + 1) -ую точку, а остальные выбросить.

Найдем закон распределения промежутка времени между соседними событиями в потоке Эрланга k -го порядка ( Эk ). Рассмотрим на оси Ot простейший поток с интервалами 1, 2,.... Величина представляет собой сумму k + 1 независимых СВ где 1, 2,..., k +1 - независимые СВ с экспоненциальным распределением.

Обозначим wk ( ) ПРВ величины для потока ( Эk ). Произведение wk ( ) d есть вероятность того, что величина примет значение между и +. Следовательно, последняя точка промежутка должна попасть на элементарный участок (, + ), а предыдущие k точек простейшего потока – получим При 0 находим точное равенство для ПРВ интервалов времени для потока Эрланга k -го порядка интервалов между событиями в потоке Эрланга Dk = ( k + 1) 2. Интенсивность порядка увеличиваются математическое ожидание и дисперсия промежутка времени между событиями, а интенсивность потока падает.

Рассмотрим, как будут изменяться характеристики потока Эрланга при увеличении k, если его интенсивность будет сохраняться постоянной.

Пронормируем величину так, чтобы ее математическое ожидание (и, следовательно, интенсивность потока) оставались неизменными. Для этого Назовем такой поток нормированным потоком Эрланга k -го порядка. ПРВ интервала между событиями Математическое ожидание mk = 1 величины, распределенной по % закону wk ( ), не зависит от k. Дисперсия интервала между событиями означает, что при неограниченном увеличении k нормированный поток Эрланга приближается к регулярному потоку с постоянными интервалами, равными 1. Это свойство потоков Эрланга дает возможность, задаваясь различными k, получать любую степень последействия: от полного отсутствия ( k = 0 ) до жесткой функциональной связи между моментами появления событий ( k = ).

5.1.4. Нормальный поток событий Поток Пальма, для которого интервалы времени распределены по нормальному закону называется нормальным потоком [21].

Следует заметить, что понятие нормального потока можно рассматривать только как приближенное. Дело в том, что при нормальном законе распределения СВ изменяется в интервале (, +), а интервал между событиями - положительная СВ, изменяющаяся в интервале (0, ). Однако интервал практически возможных значений СВ, распределенной по нормальному закону, равен (m 3, m + 3 ).

Таким образом, при условии m > 3 можно сказать, что вероятность отрицательных значений СВ практически равна нулю, так как она меньше 0, (правило «трех сигм»). Именно для этого можно использовать нормальный закон распределения. В пределе, при 0, нормальный закон распределения переходит в регулярный.

Рассмотрим связь потока Эрланга с нормальным потоком. ПРВ СВ, равной интервалу между заявками в потоке Эрланга, определяется по формуле (5.5).

Функция распределения получается интегрированием Для функции составлены специальные таблицы. При больших значениях параметра a ( a > 20 ) для вычисления R ( k, a ) можно использовать таблицы значений нормальной функции распределения:

где 0 ( x) = 5.2. Самоподобные (фрактальные) модели случайных потоков Рассмотренные пуассоновские модели потоков, когда ПРВ имеет показательный вид и входной поток обладает свойством марковости (см. главу 3), оказываются неадекватными при анализе трафика в сетях пакетной коммутации таких как Ethernet, Internet, Telnet и др. В подобных сетях На интуитивном уровне это означает, что число событий на заданном временном интервале может зависеть от числа событий, поступивших в весьма отдаленные от него интервалы времени. При этом часто процесс носит «пачечный» характер. Проведенные в последние годы исследования показали, что поведение потока данных в таких сетях связи хорошо описывается с помощью моделей, основанных на теории фракталов [25].

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

Аналогичная зависимость определяет длину некоторых рек, побережье многих островов, путь, проходимый частицей при броуновском движении, и многое другое. В качестве наглядного примера можно привести так называемый «остров Коха». На рис. 5.5 показано, как можно построить такую фигуру [24].

На первом шаге берем обычный равносторонний треугольник (рис. 5.5).

Потом на каждой стороне достраиваем по треугольнику, сторона которого в три, а значит, площадь в девять раз меньше, чем у исходного треугольника и так далее. То, что получится после бесконечного количества таких шагов, называется островом Коха. При этом длина его побережья бесконечна, поскольку на втором шаге периметр фигуры увеличится в 4/3 раза, на третьем еще в 4/3 и т. д. Это происходит потому, что каждый отрезок мы заменили ломаной, длина которой в 4/3 раза больше. Таким образом, периметр данной фигуры p = lim ( 4 3) =. При этом с помощью формул геометрической прогрессии можно убедиться, что площадь острова Коха конечна.

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

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

В крупном масштабе видим большой зубец и несколько маленьких. Увеличим маленький зубчик. То есть, по существу, увеличим маленький прямоугольник до размеров первоначального. Опять выделим маленький прямоугольник, опять увеличим и опять увидим то же самое и так до бесконечности. Это свойство, выглядеть в любом, сколь угодно мелком масштабе примерно одинаково, называется масштабной инвариантностью, а множества, которые им обладают - фракталами [24]. Можно спросить, как же характеризовать фракталы, если размеры становятся какими-то зыбкими, ненадежными и начинают зависеть от размеров линейки?

Рис. 5.6. Пример фрактальной масштабной инвариантности На это математики могут ответить просто и остроумно: «Важна не сама длина, а то, как она зависит от размеров линейки, т. е. важно некое число, называемое фрактальной размерностью». Для отрезка размерность равна 1, для квадрата - 2, для куба - 3. Для фракталов размерность - дробное число. Отсюда и само название «фрактал», происходящее от английского «fractal» - дробный, неполный, частичный. Например, для острова Коха оно лежит между 1 и 2 полоса в двумерном пространстве, т. е. уже не обычная кривая, но еще не плоскость (дробная размерность) [24].

5.2.2. Самоподобные (фрактальные) случайные процессы Представленный пример фрактала (кривая Коха) относится к классу детерминированных фракталов, т. е. когда объект непосредственно составляется из своих малых копий. В теории телетрафика для описания поведения величины нагрузки в сетях связи с пакетной коммутацией применяется класс случайных (стохастических) фракталов. В этом случае свойство самоподобия (масштабной инвариантности) наблюдается лишь «в среднем», т. е. подобными являются не сами отсчеты сигнала, а, например, его КФ или ПРВ в разных временных масштабах. Три характерные особенности самоподобных процессов выражены в медленном убывании дисперсии, долгосрочной зависимости и флуктуационном характере спектра мощности таких процессов [25].

Рассмотрим дискретную случайную последовательность отсчетов:

где xi - СВ с заданным законом распределения. Будем предполагать, что все рассматриваемые СП имеют ограниченную ковариацию B ( xi, xi + ) <, и следовательно дисперсию xi = B ( xi, xi ) <. СП будет обладать свойством самоподобия, если агрегированный процесс m -го порядка будет иметь КФ r ( m ) (k ) совпадающую с КФ r (k ) исходного СП для любых m.

При выполнении данного условия можно утверждать, что дисперсия агрегированного процесса X (m ) убывает согласно выражению т. е. дисперсия агрегированных процессов – средних выборок – уменьшается медленнее, чем величина, обратная размеру выборки. В результате в самоподобных процессах имеет место явление долгосрочной зависимости, которое приводит к расходимости КФ процесса:

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

Важнейшим параметром, характеризующим «степень» самоподобия СП, является параметр Хёрста.

Для выборочного случайного набора X j ( j = 1,..., N ) можно определить выборочное среднее выборочную дисперсию и интегральное отклонение Определим изменчивость СП на интервале N как неубывающую функцию длины интервала Хёрстом было показано, что для большинства естественных процессов при больших значениях N выполняется соотношение или иначе Величина H получила название параметра Хёрста и лежит в интервале 0.5 < H < 1.0. Для процессов, не обладающих свойством самоподобия, H = 0.5.

Для самоподобных процессов с долгосрочной зависимостью этот параметр изменяется в пределах 0.7…0.9. Параметр, который был введен выше для задания асимптотических свойств характеристик самоподобных СП, может быть выражен через параметр Хёрста:

Степень самоподобия процесса можно оценить другим способом: параметр Херста можно определить путем построения графика отношения log ( RN / S N ) в зависимости от log ( N / 2 ) при разных N и вычислить величину H как тангенс угла наклона полученной линии. Следует заметить, что полученное множество точек не будут лежать на одной линии, поэтому их следует аппроксимировать линией, например, по методу наименьших квадратов. Данная методика определения параметра Херста получила название R/S-метод [25].

R/S-метод дает лишь приближенное значение показателя Херста, поэтому для его вычисления целесообразно пользоваться несколькими методиками и сравнения полученных результатов. Рассмотрим метод определения величины H на основе периодограммного анализа.