Содержание к диссертации
Введение
Глава 1. Моделирование некоторых классов кусочно-постоянных негауссовских процессов и полей 18
1.1 Кусочно-постоянное восполнение стационарных в широком смысле процессов с узлов регулярной сетки в произвольную точку временного интервала 18
1.2. Восполнение однородного в широком смысле поля с узлов регулярной сетки в произвольную точку области 28
1.3 Восполнение однородного и изотропного ПОЛЯ 34
1.4 Кусочно-постоянное восполнение стационарных в широком смысле процессов с узлов стохастической сетки в произвольную точку временного интервала 46
Глава 2. Моделирование некоторых классов гауссовских процессов и полей с корреляционными матрицами теплицева и блочно-теплицева вида 57
2.1 Рекурсивные алгоритмы моделирования реализаций условно распределенных стационарных гауссовских процессов 57
2.2 Некоторые приемы регуляризации алгоритмов моделирования векторных гауссовских последовательностей 67
Глава 3. Численное моделирование стохастической структуры слоисто-кучевых облаков 76
3.1 Каскадная модель 77
3.2 Модели на основе метода нелинейного преобразования гауссовского поля. Спектральные модели 83
3.3. Схема скользящего суммирования. Сравнительный анализ моделей SS
Заключение 92
Литература
- Восполнение однородного в широком смысле поля с узлов регулярной сетки в произвольную точку области
- Кусочно-постоянное восполнение стационарных в широком смысле процессов с узлов стохастической сетки в произвольную точку временного интервала
- Некоторые приемы регуляризации алгоритмов моделирования векторных гауссовских последовательностей
- Модели на основе метода нелинейного преобразования гауссовского поля. Спектральные модели
Введение к работе
В связи с бурным развитием вычислительных средств
появляется возможность решать более сложные и трудоемкие
научные и прикладные задачи в различных областях науки и
техники с привлечением методов статистического моделирования,
требующие соответствующих алгоритмов для их реализации. При
решении современных задач, связанных с моделированием и
анализом атмосферных процессов и климата [14,26,31-34,35,40-
43,45,52,54,56], океанологических и гидрологических процессов [4,
22, 25, 63-65, 69], процессов переноса примесей в атмосфере [62],
при решении различного типа природоохранных задач [21], задач по
исследованию закономерностей переноса излучения в атмосфере
[30,37,75,84,85,97], в задачах статистической турбулентности
[94,95], при объединении гидротермодинамических и
вероятностных подходов к описанию реальных процессов [53,71,73,74,92], при исследовании климатических воздействий атмосферных процессов на различные строительные объекты и сооружения [13], при решении агрофизических задач [29] и т.д. используются многомерные стохастические модели атмосферных процессов. Эти задачи определяют новые требования к численным вероятностным моделям - увеличение размерности, привлечение большого объема фактической информации, учет в моделях большого числа статистических параметров. При этом для решения перечисленных выше задач необходимы также реалистичные стохастические модели атмосферных процессов, в достаточной степени адекватно описывающие характерные особенности реальных данных. В первую очередь в них должны быть учтены доступные (в рамках имеющейся информации) особенности
распределений и корреляционных связей многомерного процесса. С
другой стороны, при использовании вероятностного подхода к
решению этих задач необходимо моделировать очень большое число
реализаций процессов и полей для получения надежных оценок
рассчитываемых функционалов. В связи с этим разработка
универсальных и экономичных методов и алгоритмов
моделирования многомерных случайных процессов и полей с
заданной корреляционной структурой (в том числе с учетом
нестационарности или неоднородности) и алгоритмов
стохастического восполнения, используемых при построении стохастических моделей атмосферных процессов, является актуальной теоретической и важной практической задачей.
Основными характеристиками, используемыми при
построении численных алгоритмов моделирования случайных
процессов и полей, являются одномерные распределения,
корреляционные функции либо спектральные плотности
соответствующих процессов и полей. Для учета одномерных
распределений используются различные подходы, например
функциональные нелинейные преобразования гауссовских и негауссовских процессов и полей, моделирование кусочно-постоянных процессов и полей с использованием различных типов точечных потоков, линейные преобразования негауссовских процессов.
Чаще всего в основе методов моделирования негауссовских процессов и полей лежат методы моделирования гауссовских процессов и полей. Это главным образом связано с тем, что для гауссовских процессов и полей класс корреляционных функций (или корреляционных матриц в случае поля на ограниченной сетке) наиболее широк, он определяется неотрицательной
определенностью этих корреляционных функций или матриц. Методы моделирования гауссовских векторов с нулевым средним и заданной ковариационной матрицей R основаны на использовании линейного преобразования [7, 28, 70]
І~Аф (1)
гауссовского вектора ф = (<рх,...,<рп)т с нулевым средним, единичными дисперсиями и независимыми компонентами, где А -нижняя треугольная матрица такая, что АЛТ = R. Методы разложения симметричной положительно определенной матрицы R на произведение транспонированных друг относительно друга треугольных матриц хорошо известны (метод Холецкого) [9,79], однако при больших размерностях матрицы R, а также при ее плохой обусловленности моделирование затруднительно. В том случае, когда матрица вырождена, может быть использовано
представление гауссовского вектора в виде % = РА]/2ф, где Р -матрица из собственных векторов-столбцов ковариационной матрицы R, а Л - диагональная матрица из ее собственных чисел. Для реализации этого преобразования необходимо решить полную спектральную проблему для матрицы R. Известные итерационные методы [78, 79] и современные вычислительные средства позволяют решать эти задачи с достаточной точностью для матриц большой размерности п, например, п порядка нескольких тысяч. В приложении к статистической метеорологии такого типа задачи были рассмотрены в работе [27]. Однако в ряде случаев, когда необходимо моделировать более длинные последовательности, приведенные выше методы оказываются непригодными. В этих условиях при моделировании используют различного рода приближения, например стационарное приближение. В случае
стационарных временных рядов с ограниченным числом элементов
соответствующая ковариационная матрица
Я = (^) = (^-,]) = (^) = ^,,,/,./ = 1,...,^ имеет теплицеву структуру,
которая полностью определяется своей первой строкой. Для таких матриц эффективными оказываются алгоритмы моделирования, основанные на методе условных распределений [28,70,72,88,96], основанные на преобразованиях вида [46, 90]
где
а Ь[к] и dk определяются уравнениями [39, 46 ,90]
dk = \-?kTbT[kl (2)
Вычисления Ь[к] и dk при заданной матрице Rk осуществляются рекурсивно с помощью алгоритма Дурбина [2,10,46,86,90].
Для моделирования последовательностей большой длины используются также различные модели стационарных процессов, например, модели авторегрессии, а также смешанные модели авторегрессии-скользящего среднего. При моделировании гауссовских процессов авторегрессии к - го порядка
1=^-1 + ^2+... + ^ +dq>, (3)
параметры модели однозначно определяют корреляционную структуру временного ряда. Так, например, если параметры модели определяются уравнениями (2), и при этом расширенная матрица Rk+X положительно определена, то процесс (3) является стационарным, его корреляционная функция в первых h = l,...,k точках совпадает с первой строкой ковариационной матрицы Rk, а ее значения при h > к определяются разностным уравнением [44, 90]
В ряде случаев, когда требуется знание конкретных выражений для корреляционных функций либо спектральных плотностей в виде дробно-рациональных выражений удобно пользоваться смешанной моделью авторегрессии-скользящего среднего АРСС(,/) [2,5,87]
= biLi + Ьг^-г + + ьЛ(-к + di
Рассмотренные подходы к моделированию гауссовских
векторов и процессов обобщаются на случай многомерных
процессов дискретного аргумента, а также пространственных и
пространственно-временных скалярных и векторных гауссовских
полей на регулярных и нерегулярных сетках [81,90]. В этих случаях
вместо скалярных рассматриваются блочные ковариационные
матрицы, а для однородных, а также однородных и изотропных
полей - блочно-теплицевы матрицы с блоками, обладающими
дополнительной спецификой, определяемой свойствами
соответствующих корреляционных функций. Более конкретно эти
свойства будут описаны в первой главе данной диссертации.
Другим важным классом алгоритмов моделирования стационарных
гауссовских процессов и однородных скалярных и векторных
гауссовских полей непрерывного аргумента являются алгоритмы,
основанные на спектральном представлении случайных процессов и
полей [28,50,57-59,68,77]. В отличие от алгоритмов и моделей,
основанных на линейных преобразованиях, рассмотренных выше,
которые в пределах точности датчиков гарантируют гауссовость
моделируемого процесса или поля, спектральные модели являются
приближенными и соответствующие процессы и поля при конечном
числе используемых гармоник являются приближенными
гауссовскими и становятся гауссовскими лишь асимптотически. Соответствующие теоремы о сходимости спектральных моделей приведены в работах [11,12,50,58,]. Спектральные модели находят широкое применение при решении задач, в которых необходимо знать значение процесса или поля в произвольной точке области, например, при решении задач, связанных с переносом излучения в атмосфере [90].
При решении ряда прикладных задач большое значение имеет учет нестационарности процесса или неоднородности поля. В рамках преобразования (1) для процессов и полей дискретного аргумента общая задача моделирования нестационарных гауссовских процессов и полей решается автоматически, но требует знания соответствующей ковариационной матрицы, как правило, очень большой размерности. В ряде случаев с учетом специфики решаемой задачи такую матрицу можно задавать модельно, используя некоторые функциональные соотношения для расчета ее элементов, например, использовать прямое произведение корреляционных матриц [90]. При сравнительно небольших размерностях матрицы и при достаточно большом объеме данных о процессе ковариационная матрица может быть построена с помощью соответствующих оценок по данным наблюдений. Однако в реальных условиях эта матрица может содержать большие статистические погрешности, быть плохо обусловленной либо вырожденной, поэтому для реализации преобразований типа (I) может оказаться необходимой специальная корректировка матрицы либо специальные модификации алгоритмов для расчета матрицы преобразования [3, 50].
В ряде случаев в приложениях используются алгоритмы моделирования специальных классов нестационарных процессов и неоднородных полей. В частности к ним относятся алгоритмы моделирования периодически коррелированных гауссовских и негауссовских процессов [25], а также алгоритмы моделирования частично однородных пространственных полей [98]. Для моделирования периодически коррелированных процессов используются два типа алгоритмов. К первому относятся алгоритмы,
основанные на спектральном представлении периодически коррелированных процессов в виде
где gp (t) -стационарные процессы, образующие векторный
бесконечномерный стационарный случайный процесс [25]. В
соответствующих приближенных алгоритмах могут быть
использованы спектральные модели векторных стационарных
процессов с заданными матричными ковариационными функциями.
Другой подход к моделированию периодически нестационарных
скалярных процессов, рассмотренный, например, в [66] основан на
использовании свойств блочно-теплицевых ковариационных
матриц, в которых периодичность корреляций заложена по
определению. Алгоритм моделирования периодически
нестационарного скалярного процесса
с периодом р в гауссовском случае основан на рассмотренном выше алгоритме моделирования стационарной последовательности векторов
Ц\ > Цг і Vt >
с матричной ковариационной функцией R0,RlyR2 Здесь
fji = (т](ір +1), - - -, Т](ір + р))т, р - число значений ряда внутри интервала, соответствующего периоду, / = 1,2,...- порядковый номер интервала. Смежные к элементов в последовательности 7Jltf}2f...7Jtt... образуют случайный вектор rjw = (rjl ,...,%) с
блочно-теплицевой ковариационной матрицей
R = (Д._ Л i,j = \,,..,k. При этом размерность матриц
1!
^|i-/| ~ ^1 равна p. Первая блочная строка этой матрицы
определяет первые k матриц в матричной ковариационной функции R0,RX,R21.... Для учета негауссовости могут быть использованы
различные нелинейные преобразования, периодически зависящие от времени.
При разработке негауссовских спектральных моделей также
используется модификация метода функциональных
преобразований гауссовского процесса [83], учитывающая специфику спектральных моделей. Приближенная модификация "метода обратных функций", предназначенная для моделирования временных рядов по данным реальных наблюдений, основана на нормализации реального ряда и рассмотрена в работах Г.Г. Сванидзе [69], А.С. Марченко, А.Г. Семочкина [47] . Для построения стохастических моделей широко используются также различные функциональные преобразования гауссовских процессов [28,38]. Одним из наиболее распространенных методов моделирования негауссовских процессов и полей является метод обратных функций распределения [28,55,57], впервые предложенный З.А. Пиранашвилли [55]. Суть этого метода состоит в
следующем. Пусть < =(р...,я)г- гауссовский вектор с нулевым средним и корреляционной матрицей G = (gy), i,j = 1,...,п. Для
того, чтобы построить случайный вектор fj = (??,,..., г}п)т с корреляционной матрицей R =(^), (,,/ = 1,...,и и одномерными распределениями F^x) его компонент используется преобразование каждой компоненты вектора в виде
7,=^(^,)), С4)
где Ф(лг) - функция стандартного одномерного нормального распределения. Коэффициенты корреляции ri} и gSJ связаны между собой соотношениями
MMlZ^BL, (5)
+ со
мт = \\F[xmxw;\4y))f(x,y;g,)dxdy,
— со
а /(^^5) " двумерная плотность распределения двух нормальных
величин с нулевыми средними, единичными дисперсиями и коэффициентом корреляции gg. Функция rij=f{glJ) является
монотонной, однако интервал [-1,1] изменения коэффициента корреляции gy. в зависимости от функций распределенияFt(x)и
Fj (х) отображается в более узкий интервал [а, Ь]. Для построения вектора 77 с функциями распределения Ft(x) необходимо решить систему нелинейных уравнений r^-=-f{giJ) относительно gif по заданным /v, и, в случае, если решение существует для каждого і и у, и, при этом, матрица G из коэффициентов корреляции gtj оказывается положительно определенной, то построение вектора Ц
сводится к моделированию гауссовского вектора с с корреляционной матрицей G с последующим преобразованием каждой компоненты этого вектора по формуле (4). Если при заданных /v и Fk{x) решения уравнения^. =f(gjj) не существует, то
возможны приближенные решения, в частности, в качестве gi}
выбирается наибольшее из возможных, но такое, которое гарантирует неотрицательную определенность матрицы G. Условия
совместимости одномерных распределений и корреляций приведены в [57,60]. В работе [16] приведены таблицы одномерных распределений, которые могут быть использованы для численного моделирования негауссовских процессов.
Для моделирования негауссовских процессов и полей с
произвольным одномерным распределением и произвольной
выпуклой корреляционной функцией хорошо известны методы,
основанные на различных модификациях метода "повторений" [49],
а также методы, основанные на использовании точечных потоков
Пальма [48,51] (Г.А. Михайлов). В ряде случаев могут быть
использованы также методы моделирования негауссовских
кусочно-постоянных процессов и полей с непрерывным аргументом,
основанные на объединении методов моделирования процессов и
полей дискретного аргумента с методами моделирования процессов
и полей на точечных потоках [91]. Необходимость в таких
алгоритмах определяется тем, что модели процессов и полей с
дискретным аргументом, рассмотренные выше, несмотря на
экономичные алгоритмы, основанные на использовании
специфики теплицевых корреляционных матриц, тем не менее,
ограничены ресурсами памяти ЭВМ. Наиболее существенно это
при моделировании пространственных и пространственно-
временных полей на сетках с большим числом узлов. В этих случаях
одним из возможных подходов к решению задач большой
размерности является использование методов стохастического
восполнения случайных процессов и однородных случайных полей
дискретного аргумента с узлов регулярной сетки в заданную точку
области, рассмотренных в работах [17,18,19,36,89], которые в свою
очередь основаны на модификации известного метода повторений
[49] и метода моделирования случайных процессов с
использованием стационарных точечных потоков [48,50]. Основная специфика этих алгоритмов состоит в сохранении вероятностных свойств исходного процесса (или поля на сетке), главным образом, одномерных распределений, свойств стационарности в широком смысле для процесса, однородности и изотропности для поля.
В первой главе диссертации рассматриваются вопросы, связанные с моделированием некоторых классов негауссовских случайных процессов и полей, включающих стационарные в широком смысле кусочно-постоянные процессы, однородные в широком смысле и однородные изотропные поля, конструктивное представление которых основано на использовании специальных классов точечных потоков. Полученные результаты являются частными решениями более общей задачи кусочно-постоянного стохастического восполнения пространственных полей с узлов нерегулярной сетки в произвольную точку области. Для этой задачи автором получены простые численные алгоритмы для вычисления корреляционной функции поля по заданной дискретной корреляционной функции поля в узлах сетки.
В этой главе рассматривается также алгоритм кусочно-постоянного восполнения дискретного процесса, заданного в случайные моменты времени. Рассматривается модель, в которой эти моменты времени образуют пуассоновский поток точек. Исследована корреляционная структура таких кусочно-постоянных процессов. Получены выражения для корреляционных функций при кусочно-постоянном восполнении слева и справа, проведены тестовые расчеты. Рассмотренная задача может иметь практический интерес, т.к. измерения в случайные моменты времени - ситуация, часто встречающаяся на практике.
В главе 2. рассматриваются алгоритмы моделирования условно распределенных стационарных процессов и однородных случайных полей с использованием алгоритмов, основанных на методе условных математических ожиданий. Для гауссовского случая в пункте 2 Л построены экономичные рекурсивные алгоритмы стохастической интерполяции стационарных случайных процессов с узлов регулярной сетки на регулярную сетку с меньшим шагом. При построении алгоритма используются алгоритмы метода условных математических ожиданий для теплицевых корреляционных матриц. В пункте 2.2 рассматривается алгоритм моделирования многомерных гауссовских последовательностей с блочно-теплицевыми ковариационными матрицами [90]. В ряде случаев этот алгоритм с блочно-теплицевыми ковариационными матрицами большой размерности может оказаться вычислительно неустойчивым. В связи с этим в пункте 2.2 рассматривается способ регуляризации этого численного алгоритма. Здесь основная задача состоит в том, чтобы обеспечить устойчивость вычислений и, при этом не нарушить специфики корреляционной матрицы, т.е. сохранить стационарность векторной последовательности. Предложен способ регуляризации, который сводится к тождественной записи исходного алгоритма, но, при этом, новый алгоритм содержит параметр, варьируя которым, можно добиться устойчивости вычислений. Алгоритм реализован на ЭВМ, проведены численные эксперименты, подтверждающие его эффективность.
В главе 3. Исследованы свойства трех известных алгоритмов численного моделирования случайных полей с логнормальным одномерным распределением - каскадная модель, модель, основанная на схеме скользящего суммирования и спектральная
модель. Алгоритмы предназначены для моделирования слоисто-кучевой облачности применительно к задачам переноса излучения в атмосфере. Для каскадной модели выписываются моделирующие формулы, и показывается, что в этом случае полученное поле не является однородным и дает менее реалистичную картину, по сравнению с остальными двумя алгоритмами моделирования. Все три алгоритма численно реализованы, проводится их сравнительный анализ.
Результаты, изложенные в диссертационной работе, докладывались и обсуждались на Семинаре отдела «Статистическое моделирование в физике» ИВМиМГ СО РАН, на 4-х конференциях молодых ученых ИВМиМГ СО РАН (1996-2004 гг.), на конференции IRS 2000: Current Problems in Atmospheric Radiation proceeding of the International Radiation Symposium St.Peterburg, Russia, 24-29 july 2000.
Основные результаты диссертации опубликованы в 7 печатных работах [17-20,36,89,93].
Автор выражает признательность научному руководителю д.ф.-м.н. В. А. Огородникову, члену-корреспонденту РАН Г. А. Михайлову за поддержку и постоянное внимание к работе, профессору СМ. Пригарину и к.ф.-м.н. А.В. Протасову за полезные обсуждения и сотрудничество.
!7
Восполнение однородного в широком смысле поля с узлов регулярной сетки в произвольную точку области
В работе [36] рассмотрен следующий алгоритм кусочно-постоянного восполнения негауссова однородного в широком смысле поля с узлов регулярной сетки щ = {(х0 ± hi,x0 ± hj)\ і, J = 0,1,---; xQ = 0} (для простоты будем рассматривать сетку и?, = {(х0 ± /,xQ ±j); itj = 0,1,...; х0 = 0} с целочисленными координатами узлов) в произвольную точку области D на плоскости.
Пусть (i,j)=%y однородное в широком смысле негауссово двумерное поле в узлах сетки озх с одномерным распределением F{u) и корреляционной функцией
Рассмотрим процедуру восполнения поля . с узлов сетки й ! в произвольную точку (х,у) области D (Рис, 2). (Id) Независимо выбираются две случайные величины ах и а2, равномерно распределенные в интервале [0,1]. На каждой оси х и у выбираются точки i±avj±a2, где і и j - координаты сетки на осях х и у. Перпендикулярно осям через точки і ± щ, j ± а2 проводятся прямые, разбивающие всю плоскость на квадраты со сторонами единичной длины. Выберем эти квадраты таким образом: Tif(.xty) = {xty:i-l + a1 х i + ai\j-\ + a2 y j + a2)
Утверждение 1.2.1. Одномерное распределение поля %(х,у) в произвольной точке (х,у) совпадает с F(w), и поле является однородным в широком смысле с корреляционной функцией Л(Дх,Ду) = 2 К(Ах-т,Ау-п)Птгг, (1.2.1) w=-w я=-» где О, Ддг 1и/шД » 1 Доказательство, Получим соотношение (1,2.1). Вначале рассмотрим случай, когда . не коррелированны. Для простоты будем рассматривать случай, когда Л = 0, =1. (1.2.3) Покажем, что корреляционная функция поля %(х,у)г полученного в результате преобразования (ld)-(2d) поля . имеет вид (1.2.2).
Действительно, коэффициент корреляции между точками (x0,yQ) и (х0+ Ах,у0 + Ау) можно представить в виде (здесь рассматриваем Ддг 1, Ау 1, т.к. иначе коэффициент корреляции будет равен нулю в связи с некоррелированностью %у): R(x0,yQ;xQ + &х Уа + АУ) = Щ(х0,у0){;(х0 + Ax,yQ + Ay) = = P(Atj)M +Д4+,Ж й+І + P(AMj)M4lv + P(AMJ+i)MClj+„ 2 4) где Р(Ау), P(AiJ+l), P(At+ij), P(Aj+lJ+l) - вероятности событий, состоящих в том, что точки (JT0,J 0), (x0 + Axty0 + Ay) принадлежат квадратам TfftT ltTMJ и 7]+1 +1, т.е. событий вида 4 = (іхо Уо)є Ть (хо + ДХ Уо+ АУ)е Tij) зо Пусть х 0 = Е(х0), у 0=Е(уй). Учитывая равномерность и независимость величин ах и а2, условия (1.2.3) и некоррелированность величин ff получаем, что при Ах 0, Ау О К( Уо\Хо+& Уъ+Ьу) = Р{А + Р(А х) + Р(АМ1) + Р{Ам ) = = (1-х 0 Ах)(\-у 0-Ау) + (1-х 0-Ах)у о+х 0(1-у 0-Ау) + х 0у = = 1 - Ах Ау 4- АхАу, при Дх 0, Ау О R(x0,y0;x0+Ax,y0 + Ay) = = (\-x Q-Ax)(\ y o) + (\-x Q-Ax)(y o + Ay) + x o(l-y,0) + x oy o = = 1 - Ах + Ау - ДсДу, при Ах 0, Ау О Я(х0,л; 0 + Дх,7о + АУ) = = (1- Х1- -Ау) + (1- ) ;+(л 0+АхХ1- -Ау) + К+Аї) = = 1 + Дх - Ду - ДхДу, при Дх 0, Ду 0 ACWo; х0 + Дх,у0 + Ду) = = (1- )(1- ) + (1- )( +Д ) + К + ДгХ1- ;) + К+Дл-)( +Ду) = = 1 + Ах + Ау + ДхДу, что и доказывает формулу (1.2.1), в которой в данном случае может быть только одно слагаемое при Rm l, т.к. не коррелированны.
Рассмотрим случай, когда значения поля s в узлах сетки коррелированны. Пусть і?тл, m,w = 0,±l,±2,... - корреляционная функция поля .. Коэффициент корреляции между произвольными точками (х0,у0), (х0+Ах,у0+Ау) поля (х,у) имеет вид R(x0,y0\xo+Ax,y0 + Ду) = = 11 0) 0+ 0 )6 + i+I y+1 +ІЕ%.л).К+ -«-і.л+ -«) А+!і„ + =/ /= i+l J+l +ЕІ% о) ( + Дж-т, + Ду-и -1) є Гя)йтп+1 + + 2 0.). + Дг-т-1, 0 + Ду-и-1) є Гй)Дт+1іЯ+І = =і t=j = ( Дх - m, Ду - n)Rmn -f K(Ax m 1, Ду - n)Rm+ifl + +K(Ax-m,Ay-n-V)Rmti+i+K(Ax-m-\,Ay-m-l)R/n+ln+l, где m = ( Дяг), и = E(Ay). Полученное соотношение может быть представлено в виде (1.2.1) [36], откуда непосредственно следует однородность поля в широком смысле.
Кусочно-постоянное восполнение стационарных в широком смысле процессов с узлов стохастической сетки в произвольную точку временного интервала
При каждом фиксированном разбиении полученное поле является сеточно-однородным и изотропным в смысле определения (1.3.1). При измельчении шага сетки все большее число узлов ложится в некоторой заданной окрестности є окружности, проведенной из фиксированного узла (/,/) через любой другой узел (i\f) исходной сетки. Среди корреляций между значениями поля в узле (i,j) и узлами в окрестности є равными оказываются только те, которые относятся к узлам, равноудаленным от точки (i,j) в силу алгоритма построения поля, который обеспечивает выполнение соотношения (1.3.1). В частности, коэффициент корреляции R .+mpW из (1.3.6) относится к узлам исходной сетки при т = 2N с расстоянием между ними, равным шагу сетки & ,, а Л;+/2 .+и/2 из (1.3.7) при m0 =round{2N Н2) (где round{.) - ближайшее целое число) соответствует коэффициенту корреляции между значениями поля в узлах, расстояние между которыми наиболее близко к шагу исходной сетки сог Так, в частности, при N = 10 R, n-,nn = 0-086 н-0.41Д№, + 0.5Д,,,,,.,. (1.3.8) Рассмотрим, например, исходное однородное изотропное в смысле определения (1.3.1) поле Щ на сетке щ с заданной корреляционной функцией Rmn = ехр(-а(т2+л2)). С помощью алгоритма (lf)-(2f) построим поле ,]!.N при N = \0. На Рис. 3 приведены графики зависимости коэффициентов корреляции R,... и R. ,„v . ,„ от а.
Видно, что при малых а, т.е. когда /?fflrt достаточно медленно убывает, эти корреляции достаточно близки и в этом смысле сеточное поле %}yN близко к непрерывному однородному и изотропному полю. При больших а эти корреляции существенно отличаются.
Пример. Рассмотрим еще один способ восполнения однородного и изотропного поля в смысле определения 1.3,1). Пусть, как и раньше - однородное и изотропное в смысле определения 1.3.1) двумерное негауссово поле в узлах сетки (У, = {(д:0 ± /, х0 ± j); iy j = 0,1,...; xQ = 0} с одномерным распределением F{u) и корреляционной функцией R(\i-k\t\j-l\) R{mtn) = Raet. Рассмотрим следующий способ восполнения: (lg) Для каждого /,/= 0,±1,+2,... в узлах (/ + 1/2;у +1/2) значение поля %)!}} полагаем равным одному из соседних узлов О »Л» Р»У + 1)і 0 + 1 У)5 (/ + 1,7+1) с вероятностью 1/4. (2g) Для каждого /,/ = 0,±1,±2,.„ в узлах (/,j + 1/2)значение поля Е,)!}} полагаем равным одному из соседних узлов (/,/), (/,/ + 1), (/-1/2,/ + 1/2), (/ + 1/2,7 + 1/2) с вероятностью 1/4, а в узлах (/+1/2,/)значение поля Ё } полагаем равным одному из соседних узлов(/,7), (/ + 1,/), (/ + 1/2,/-1/2), (/ + 1/2,7 + 1/2) с вероятностью 1/4. В узлах / = 0,+1,+2,..., / = 0,+1,+2,... значения поля 2 полагаем равным значениям поля в тех же узлах.
При этом нужно разыгрывать только две независимые случайные величины, равномерно распределенные на интервале [0,1]: одну для шага (lg), другую для шага (2g).
Утверждение 1.3.4. Поле .2, полученное в результате процедуры восполнения (lg)-(2g) является однородным и изотропным в смысле определения (1.3.1). До ка зател ь ств о.
По построению, поле, полученное после шага (lg), совпадает по полученной сетке с полем, построенным в результате процедуры восполнения (lf)-(2f), значит по доказанному утверждению 1.3.2 после шага (lg) поле однородно и изотропно в смысле определения (1.3.1). А шаг (2g) соответствует шагу (lg) для сетки coh ={(x0±hn,xQ±hm);n,m = OX--- ,Xo =} с шагом сетки h = V2/2, повернутой на угол, равный я74, следовательно, опять получаем однородное и изотропное поле в смысле определения (1.3.1). Утверждение доказано.
Отметим, что процедура восполнения (lg)-(2g), в отличие от процедуры (lf)-(2f), уже не эквивалентна процедуре восполнения (ld)-(2d) однородного в широком смысле поля с узлов целочисленной сетки по значениям однородного изотропного в смысле определения (1.3.1) поля с корреляционной функцией
Rmn в точку области (/ + к/2, j + //2), /, j = 0, ± 1, ±2,..., к, I = 0,1. Корреляционная функция поля, полученного в результате процедуры восполнения (lg)-(2g) имеет более сложную структуру и R+1 к + і N зависит уже не от четырех значений: Rij RM,/ Ru+i RMj+i как ПРИ проЦеДУРе восполнения (lf)-(2f), а от шестнадцати Rgb: a = i-1,/,/ + 1,/+2, b j - \,j,j + \,j + 2.
Таким образом, в данном пункте предложены и исследованы некоторые возможные способы стохастического восполнения однородных и изотропных в смысле определения (1.3.1) негауссовых полей, заданных на регулярной сетке с шагом h = 1, на регулярную сетку с шагом l/2jV,iV l с сохранением свойств однородности и изотропности в смысле определения (1.3.1).
Некоторые приемы регуляризации алгоритмов моделирования векторных гауссовских последовательностей
Процедура моделирования условных реализаций вектора (2.1.16) с использованием соотношения (2.1.12) принципиально не отличается от предыдущей, поскольку построение вектора zn осуществляется с помощью тех же формул (2.1.14) и при тех же самых размерностях векторов хп,у„, только векторы и fj моделируются с корреляционной матрицей, которая в данном случае определяется строкой (2.1.18), а их разбиение на подвекторы производится в соответствии с (2.1.17). Основное различие состоит в последней операции умножения вектора 1п на матрицу R12, Матрица Rn в данном случае - блочная и имеет вид R\2 (2.1.19) где ЯДя] - блоки размерности (р — 1)х«, которые можно представить
Таким образом, полученная процедура для моделирования условно распределенных стационарных гауссовских последовательностей не содержит разложения матрицы Rn2 на произведение двух треугольных (на практике эта матрица может иметь достаточно большую размерность, обозначим ее через N). Отметим, что для реализации процедуры разложения необходимо хранить в памяти ЭВМ порядка N2 /2 чисел. Вместо операции разложения получена рекурсивная процедура, в которой необходимо хранить в памяти порядка 3N чисел. Для реализации формулы (2.1.20) необходимо хранить в памяти первую строку совместной корреляционной матрицы вектора (2.1.16).
Рассмотренные алгоритмы по существу являются алгоритмами моделирования ошибок оптимальной интерполяции гауссовских стационарных процессов с узлов редкой регулярной сетки в узлы более густой сетки. Обобщение на случай негауссовских процессов может быть сделано на основе метода «обратных функций» [28,55].
Некоторые приемы регуляризации алгоритмов моделирования векторных гауссовских последовательностей.
Рассмотрим алгоритм моделирования стационарно связанных гауссовских р -мерных векторов , ,..., с заданной вещественной ковариационной матрицей блочно-стационарного (или блочно-теплицева) вида R ( ) Rl л-1 и-2 R (2.2.1) Ді-l Rn-2 где Rk,k = 0,...,n-\ Rn матрицы pxp- Вещественная положительно определенная матрица R,n, - симметрична. Но ее блочные элементы могут быть несимметричными. Рассмотрим блочную ковариационную матрицу R )f связанную с матрицей Я(/)) вида (2.2.1) соотношением где -блочная матрица перестановок, а I - единичная матрица размерности РхР- В общем случае несимметричные блоки Rk k = 0,...tn — i матрицы R связаны с несимметричными блоками Rk матрицы R,n, соотношением Rk — Rk (при этом RW = RW R(") = " о = » о = Таким образом, матрица /?(rt) имеет следующий вид В случае если матрицы Rk несимметричны, матрица R n) играет вспомогательную роль при построении алгоритма моделирования. Рассмотрим известный алгоритм моделирования векторов ,х,Е,г,...,Е,п методом условных математических ожиданий [90]. В соответствии с этим методом векторы ( , = 1,.,.,72 вычисляются рекурсивно по схеме где фІУф2,— фп- независимые гауссовские векторы размерности р, такие, что Мфкф[ = Ір) Мфкф]" -,к 1. Здесь 0 - нулевая матрица размерности рхр, h(k) = 0=1 к ) В[к] = (В([к], Втк[к])т, = 1,...,л-1, В;[к] - матрицы размера рхр, а С(. - нижние треугольные матрицы размера рхр. Численная реализация схемы (2,2.4) требует вычисления на каждом шаге матричного вектора В\1с\ и остаточной ковариационной матрицы Qk = СкСк . Одним из возможных методов вычисления В[к] для к— {,...,п \ является решение уравнений вида RwB[k] = Rk, (2.2.5) а матрица Qk связана с /?(i) и В[к] следующим соотношением: Qk=RQ-BT[k]Rik)B[k].
Модели на основе метода нелинейного преобразования гауссовского поля. Спектральные модели
Метод нелинейного преобразования гауссовского поля был предложен в работе [55]. Применим его для моделирования однородного поля с логнормальным одномерным распределением.
Однородное случайное поле т(х,у) с логнормальным одномерным распределением с параметрами с, и и ковариациями r{t,s) можно моделировать (методом обратной функции распределения) по формуле r(x,y) = exp( Tw(x,y) + c), (3.2.1) где w(x,j )- однородное случайное гауссовское поле с нулевым средним, единичной дисперсией и корреляционной функцией p(t,s): p(t,s) = R-Fl(r(t,s)), +00 +00 RF(p) = j / (0( )) (0(7)) ( 7, Г +Пг-2р 2лгл/Г р ехр ФАїл) = V 2(1- ) где F - функция логнормального распределения с параметрами с, а, а Ф - функция стандартного нормального распределения, (Ф(йНхрК+с). При этом ковариации гауссовского и логнормального полей связаны следующим образом: ехр( г2р)-1 RF\P)- 1 2N -. explcr 1-І ;1(r)= 9=or 2ln[1+r(exp(0 2)"1)] где RF(P)- 2 mF и rF - математическое ожидание и дисперсия логнормального распределения с параметрами с, ег: mF - ехр — + с 2 сг = (ехр(с72)-і)ехр(ст2 + 2с). V -6 У При о 1 искажение корреляции незначительно: max рє[0,і] RF{P)-P олз, и для простоты мы будем им пренебрегать, полагая p(t,s)-r(tts).
Т.о. приближенное моделирование случайного поля тух, у) с логнормальным распределением и степенным спектром можно свести к моделированию гауссовского поля w(x,y) со степенным спектром и затем к преобразованию (3.2.1). Для моделирования гауссовского поля воспользуемся спектральным методом или методом скользящего суммирования. Спектральные модели Приведем основные моделирующие формулы, которые использовались [57, 61] для численного моделирования однородного гауссовского случайного поля с нулевым математическим ожиданием и спектрами вида ЛЛЛ) [л/Л2+ ] -5/3 В качестве численной аппроксимации использовалось выражение w4 ) = E { C0KV+v )+ sin(V+v )} и і = 1,..., п, j = -n,...-\,1,...,n, где .., тх,. - независимые стандартные нормальные случайные величины. Это выражение можно записать в эквивалентном виде, более удобном для моделирования / = 1,...,И, j =-/7,...-1,1,..., П, здесь СГу, Д.. - независимые случайные величины, равномерно распределенные на интервале [ОДJ. Для спектральной плотности f(X,v) = c0\Xv{5 \ Хє[А,В], ve[45]UbM], где 0 А В, с0 —const, интервал [ ,5] разбивался на п под интервалов А-Ъ ЪХ „. Ьп =В, величины Су определялись из равенств С-j — C- _jy I — 1,..., r?, J ij-» ) "j S5 а величины /L, v.. выбирались следующим образом: vv = (bH +bj)/2, /, у = 1,..., и, =-( ,.,+ )/2, / = l,...,w у = -!,...,-«, (нерандомизированная модель), или моделировались случайным образом в соответствующих интервалах: vy e[bj_l9bj), /,У = 1,...,я, v.. є[- ,- ), г = 1,...,/і, у = -1,...,-и, независимо друг от друга по плотности пропорциональной \Л\ (рандомизированная модель). При этом для моделирования случайных величин /L можно воспользоваться формулой =[пл-4+а-п) ]"1" . =-5/з, где у.. - случайные величины, равномерно распределенные на интервале [ОД]. Аналогично моделируются величины у.(..
Нерандомизированная модель является гауссовской, но воспроизводит корреляционную функцию приближенно, т.к. имеет дискретный спектр. Рандомизированная модель передает корреляционную функцию и спектр точно, но является лишь асимптотически гауссовской. Рандомизированную модель целесообразно использовать, когда независимые реализации случайных полей моделируются многократно.
Рассмотрим теперь однородное изотропное гауссовское поле w(x,y) на плоскости с нулевым средним, корреляционной функцией и о (J0 - функция Бесселя первого рода) и радиальной спектральной плотностью g(r)=c0r-5f\ гє[А,В], 0 А В. Разобьем интервал A=b0 Ь{ .., Ьп— В и в качестве модели изотропного поля w рассмотрим [57] /7=і и=і (3.2.2) xcos[(xrncosam+yr„sme wa) + 2xpmn]9 где К cl=cQ\g{r)dr, а}пт=к{т-уп)1Мп, ccnm,j3„m,Yn -независимые случайные величины, равномерно распределенные на [0,1]. Спектральная модель (3.2.2) состоит из ,Мп гармоник и соответствует разбиению спектрального пространства на N колец и п -го кольца на Мп секторов. Для нерандомизированной спектральной модели величины гп выбираются детерминированным образом в интервалах \рп_{,ЬЛ (например, берутся центры этих интервалов), а для рандомизированной спектральной модели величины гп моделируются случайным образом на интервалах [Ьп_:,Ьп с плотностью вероятности g(r)Jcn 3.3. Схема скользящего суммирования. Сравнительный анализ моделей. Схема скользящего суммирования удобна для моделирования гауссовских случайных полей с дискретным аргументом [76,28]. Приведем моделирующие формулы для однородного гауссовского случайного поля wnm на плоскости с нулевым средним и единичной дисперсией: W пт = Z А р«-к -р 4 со к я —К —К 1 w ft »v) = -—г Кпт cos(A« + vm) = I {2пУ AkJcos(Zk + vj) k,J=\ a(Z,v) = 4.cos(/U + vy) = fan)2 f {X,v),_ -і л к AkJ= J" J Jflr( ,v)cos(A& + v/)fikWv. 9 (2 ) -Л1 -Я Здесь AkJ - коэффициенты скользящего суммирования, р -независимые гауссовские случайные величины с нулевым средним и единичной дисперсией, Ктя- корреляционная функция поля wnm, fA(X,v), X,v є[-Л",Л") - спектральная плотность случайного поля и a(A,v) - частотная характеристика преобразования скользящего суммирования.
