Численные алгоритмы моделирования и стохастического восполнения случайных процессов и полей

В связи с бурным развитием вычислительных средств появляется возможность решать более сложные и трудоемкие научные и прикладные задачи в различных областях науки и техники с привлечением методов статистического моделирования, требующие соответствующих алгоритмов для их реализации. При решении современных задач, связанных с моделированием и анализом атмосферных процессов и климата [14,26,31−34,35,4043,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] и т. д. используются многомерные стохастические модели атмосферных процессов. Эти задачи определяют новые требования к численным вероятностным моделям — увеличение размерности, привлечение большого объема фактической информации, учет в моделях большого числа статистических параметров. При этом для решения перечисленных выше задач необходимы также реалистичные стохастические модели атмосферных процессов, в достаточной степени адекватно описывающие характерные особенности реальных данных. В первую очередь в них должны быть учтены доступные (в рамках имеющейся информации) особенности распределений и корреляционных связей многомерного процесса. С другой стороны, при использовании вероятностного подхода к решению этих задач необходимо моделировать очень большое число реализаций процессов и полей для получения надежных оценок рассчитываемых функционалов. В связи с этим разработка универсальных и экономичных методов и алгоритмов моделирования многомерных случайных процессов и полей с заданной корреляционной структурой (в том числе с учетом нестационарности или неоднородности) и алгоритмов стохастического восполнения, используемых при построении стохастических моделей атмосферных процессов, является актуальной теоретической и важной практической задачей.

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

Чаще всего в основе методов моделирования негауссовских процессов и полей лежат методы моделирования гауссовских процессов и полей. Это главным образом связано с тем, что для гауссовских процессов и полей класс корреляционных функций (или корреляционных матриц в случае поля на ограниченной сетке) наиболее широк, он определяется неотрицательной определенностью этих корреляционных функций или матриц. Методы моделирования гауссовских векторов с нулевым средним и заданной ковариационной матрицей Я основаны на использовании линейного преобразования [7, 28, 70].

1 = Ау (1) гауссовского вектора ф = (<�р1,.,<�рп)т с нулевым средним, единичными дисперсиями и независимыми компонентами, где Анижняя треугольная матрица такая, что ААТ = Я. Методы разложения симметричной положительно определенной матрицы Я на произведение транспонированных друг относительно друга треугольных матриц хорошо известны (метод Холецкого) [9,79], однако при больших размерностях матрицы Я, а также при ее плохой обусловленности моделирование затруднительно. В том случае, когда матрица вырождена, может быть использовано представление гауссовского вектора в виде %=РА½ф, где Рматрица из собственных векторов-столбцов ковариационной матрицы Я, а, А — диагональная матрица из ее собственных чисел. Для реализации этого преобразования необходимо решить полную спектральную проблему для матрицы Я. Известные итерационные методы [78, 79] и современные вычислительные средства позволяют решать эти задачи с достаточной точностью для матриц большой размерности п, например, п порядка нескольких тысяч. В приложении к статистической метеорологии такого типа задачи были рассмотрены в работе [27]. Однако в ряде случаев, когда необходимо моделировать более длинные последовательности, приведенные выше методы оказываются непригодными. В этих условиях при моделировании используют различного рода приближения, например стационарное приближение. В случае стационарных временных рядов с ограниченным числом элементов соответствующая ковариационная матрица.

Л = = = = = имеет теплицеву структуру, которая полностью определяется своей первой строкой. Для таких матриц эффективными оказываются алгоритмы моделирования, основанные на методе условных распределений [28,70,72,88,96], основанные на преобразованиях вида [46, 90] вп1 = о, д, где.

1 0 0 1 0 • •• 0.

1 0 0 •• 0.

— Ь1п-1] 1 0 0 •• • ¿-п, а Ь[к] и с1к определяются уравнениями [39, 46, 90].

ЯкЬ[к] = гл* с1к=-?{Ът[к}, (2).

Ь[к] = (Ь1[к],., Ькт) Т, гк=(г13., гк) т.

Вычисления Ь[к] и с1к при заданной матрице Кк осуществляются рекурсивно с помощью алгоритма Дурбина [2,10,46,86,90].

Для моделирования последовательностей большой длины используются также различные модели стационарных процессов, например, модели авторегрессии, а также смешанные модели авторегрессии-скользящего среднего. При моделировании гауссовских процессов авторегрессии к — го порядка =, + Ъ£&bdquo-2 +. + Ьк4к +с!(Р[ (3) параметры модели однозначно определяют корреляционную структуру временного ряда. Так, например, если параметры модели определяются уравнениями (2), и при этом расширенная матрица ЯА+1 положительно определена, то процесс (3) является стационарным, его корреляционная функция в первых к = ,., к точках совпадает с первой строкой ковариационной матрицы, а ее значения при к > к определяются разностным уравнением [44, 90] ги=ЬхГк+Ьггьх+. + Ькгкк.

В ряде случаев, когда требуется знание конкретных выражений для корреляционных функций либо спектральных плотностей в виде дробно-рациональных выражений удобно пользоваться смешанной моделью авторегрессии-скользящего среднего АРСС (&-,/) [2,5,87] = + Ь2^-2 + • • •+ък%, к + с1х (р,+а2(р (х +.+.

В работе [64] приведены таблицы выражений для корреляционных функций АРСС (2,1), которые удобно использовать при моделировании гауссовских стационарных последовательностей. Дробно-рациональные спектральные плотности часто [46,90] используют для аппроксимации спектральных плотностей, оцененных по реальным рядам. Для моделирования однородных гауссовских двумерных полей дискретного аргумента часто используется алгоритм, основанный на двумерной модели скользящего среднего [76], а также модели полей марковского типа на регулярной сетке [23]. В ряде случаев для моделирования гауссовских однородных и изотропных полей с корреляционными функциями специального вида используют различного типа стохастические дифференциальные уравнения в частных производных [15].

Рассмотренные подходы к моделированию гауссовских векторов и процессов обобщаются на случай многомерных процессов дискретного аргумента, а также пространственных и пространственно-временных скалярных и векторных гауссовских полей на регулярных и нерегулярных сетках [81,90]. В этих случаях вместо скалярных рассматриваются блочные ковариационные матрицы, а для однородных, а также однородных и изотропных полей — блочно-теплицевы матрицы с блоками, обладающими дополнительной спецификой, определяемой свойствами соответствующих корреляционных функций. Более конкретно эти свойства будут описаны в первой главе данной диссертации. Другим важным классом алгоритмов моделирования стационарных гауссовских процессов и однородных скалярных и векторных гауссовских полей непрерывного аргумента являются алгоритмы, основанные на спектральном представлении случайных процессов и полей [28,50,57−59,68,77]. В отличие от алгоритмов и моделей, основанных на линейных преобразованиях, рассмотренных выше, которые в пределах точности датчиков гарантируют гауссовость моделируемого процесса или поля, спектральные модели являются приближенными и соответствующие процессы и поля при конечном числе используемых гармоник являются приближенными гауссовскими и становятся гауссовскими лишь асимптотически. Соответствующие теоремы о сходимости спектральных моделей приведены в работах [11,12,50,58,]. Спектральные модели находят широкое применение при решении задач, в которых необходимо знать значение процесса или поля в произвольной точке области, например, при решении задач, связанных с переносом излучения в атмосфере [90].

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

В ряде случаев в приложениях используются алгоритмы моделирования специальных классов нестационарных процессов и неоднородных полей. В частности к ним относятся алгоритмы моделирования периодически коррелированных гауссовских и негауссовских процессов [25], а также алгоритмы моделирования частично однородных пространственных полей [98]. Для моделирования периодически коррелированных процессов используются два типа алгоритмов. К первому относятся алгоритмы, основанные на спектральном представлении периодически коррелированных процессов в виде.

00 2 Я" Е £&bdquo-(')ехр (|— р=-СО ^ где (г)-стационарные процессы, образующие векторный бесконечномерный стационарный случайный процесс [25]. В соответствующих приближенных алгоритмах могут быть использованы спектральные модели векторных стационарных процессов с заданными матричными ковариационными функциями. Другой подход к моделированию периодически нестационарных скалярных процессов, рассмотренный, например, в [66] основан на использовании свойств блочно-теплицевых ковариационных матриц, в которых периодичность корреляций заложена по определению. Алгоритм моделирования периодически нестационарного скалярного процесса с периодом р в гауссовском случае основан на рассмотренном выше алгоритме моделирования стационарной последовательности векторов.

ЧАг" —— с матричной ковариационной функцией Здесь.

77. = (7](1р + 1),., Г}(1р +р))т, рчисло значений ряда внутри интервала, соответствующего периоду, / = 1,2,.- порядковый номер интервала. Смежные к элементов в последовательности 77, 772,.77,.,. образуют случайный вектор = с блочно-теплицевой ковариационной матрицей = (&. А /, у = 1 При этом размерность матриц Rj равна p. Первая блочная строка этой матрицы определяет первые к матриц в матричной ковариационной функции RQ, Rl, R2, — Для учета негауссовости могут быть использованы различные нелинейные преобразования, периодически зависящие от времени.

При разработке негауссовских спектральных моделей также используется модификация метода функциональных преобразований гауссовского процесса [83], учитывающая специфику спектральных моделей. Приближенная модификация «метода обратных функций», предназначенная для моделирования временных рядов по данным реальных наблюдений, основана на нормализации реального ряда и рассмотрена в работах Г. Г. Сванидзе [69], A.C. Марченко, А. Г. Семочкина [47]. Для построения стохастических моделей широко используются также различные функциональные преобразования гауссовских процессов [28,38]. Одним из наиболее распространенных методов моделирования негауссовских процессов и полей является метод обратных функций распределения [28,55,57], впервые предложенный З. А. Пиранашвилли [55]. Суть этого метода состоит в следующем. Пусть? = ?л)ггауссовский вектор с нулевым средним и корреляционной матрицей G-{gij), i, j = ,., n. Для того, чтобы построить случайный вектор rj = {rjx,., rjn) T с корреляционной матрицей R = {rij), i, j = ,., п и одномерными распределениями F^x) его компонент используется преобразование каждой компоненты вектора? в виде щ = (4) где Ф (дг) — функция стандартного одномерного нормального распределения. Коэффициенты корреляции г^ и giJ связаны между собой соотношениями где 00 00 а Их>Уёу) ~ двумерная плотность распределения двух нормальных величин с нулевыми средними, единичными дисперсиями и коэффициентом корреляции gij. Функция гу является монотонной, однако интервал [-1,1] изменения коэффициента корреляции gij в зависимости от функций распределения/^(х)и.

FJ (x) отображается в более узкий интервал Для построения вектора г с функциями распределения /*](*) необходимо решить систему нелинейных уравнений относительно по заданным г<�р и, в случае, если решение существует для каждого / и у, и, при этом, матрица О из коэффициентов корреляции оказывается положительно определенной, то построение вектора г сводится к моделированию гауссовского вектора с с корреляционной матрицей О с последующим преобразованием каждой компоненты этого вектора по формуле (4). Если при заданных и ^(х) решения уравнения^. = не существует, то возможны приближенные решения, в частности, в качестве выбирается наибольшее из возможных, но такое, которое гарантирует неотрицательную определенность матрицы С7. Условия совместимости одномерных распределений и корреляций приведены в [57,60]. В работе [16] приведены таблицы одномерных распределений, которые могут быть использованы для численного моделирования негауссовских процессов.

Для моделирования негауссовских процессов и полей с произвольным одномерным распределением и произвольной выпуклой корреляционной функцией хорошо известны методы, основанные на различных модификациях метода «повторений» [49], а также методы, основанные на использовании точечных потоков Пальма [48,51] (Г.А. Михайлов). В ряде случаев могут быть использованы также методы моделирования негауссовских кусочно-постоянных процессов и полей с непрерывным аргументом, основанные на объединении методов моделирования процессов и полей дискретного аргумента с методами моделирования процессов и полей на точечных потоках [91]. Необходимость в таких алгоритмах определяется тем, что модели процессов и полей с дискретным аргументом, рассмотренные выше, несмотря на экономичные алгоритмы, основанные на использовании специфики теплицевых корреляционных матриц, тем не менее, ограничены ресурсами памяти ЭВМ. Наиболее существенно это при моделировании пространственных и пространственно-временных полей на сетках с большим числом узлов. В этих случаях одним из возможных подходов к решению задач большой размерности является использование методов стохастического восполнения случайных процессов и однородных случайных полей дискретного аргумента с узлов регулярной сетки в заданную точку области, рассмотренных в работах [17,18,19,36,89], которые в свою очередь основаны на модификации известного метода повторений [49] и метода моделирования случайных процессов с использованием стационарных точечных потоков [48,50]. Основная специфика этих алгоритмов состоит в сохранении вероятностных свойств исходного процесса (или поля на сетке), главным образом, одномерных распределений, свойств стационарности в широком смысле для процесса, однородности и изотропности для поля.

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

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

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

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

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

Результаты, изложенные в диссертационной работе, докладывались и обсуждались на Семинаре отдела «Статистическое моделирование в физике» ИВМиМГ СО РАН, на 4-х конференциях молодых ученых ИВМиМГ СО РАН (1996;2004 гг.), на конференции 1RS 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].

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

Заключение

.

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

2. Построен алгоритм восполнения однородного изотропного поля, заданного на регулярной сетке на сетку с меньшим шагом, сохраняющий однородность и изотропность.

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

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

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

6. Исследованы и численно реализованы известные алгоритмы моделирования слоисто кучевой облачности: каскадная модель, модель, основанная на схеме скользящего суммирования и спектральная модель. Доказано, что в каскадной модели полученное поле не является однородным. Показано, что каскадная модель дает менее реалистичную картину, но является менее трудоемкой по сравнению с двумя другими моделями.