Численное моделирование сильных и длительных понижений температуры Текст научной статьи по специальности «Математика»

Численное моделирование сильных

__ «J» _ *

и длительных понижении температуры*

К. В. Деренок Новосибирский государственный университет, Россия e-mail: k.derenok@cft.ru

In this paper, we consider the algorithms for numerical modeling of periodically correlated random discrete time processes. As an application of the considered algorithms, a numerical study of statistical characteristics of strong and long-lasting temperature decreases is investigated.

Введение

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

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

* Работа выполнена при финансовой поддержке программы Президента “Ведущие научные школы”

© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2008.

(1)

(грант № НШ-4774.2006.1).

где = (С/,1, С/,2, • • • , Й,р)Т5 г — номер суток в месяце, ] — номер года. Поскольку измерения проводились четыре раза в сутки в равноотстоящие моменты времени в интервале времени месячной длины на протяжении 47 лет, то, соответственно, р = 4 — число измерений в сутки, п = 31 — число дней в месяце, т = 47 — число лет наблюдений.

Модель строится в гауссовском приближении. Такое приближение часто используется для описания реального характера временных рядов приземной температуры воздуха [3]. Поскольку температура воздуха существенно зависит от времени суток, суточный ход температуры воздуха может принципиально повлиять на длину периодов понижения температуры. В связи с этим в работе используется модель периодически коррелированных процессов, что позволяет более точно рассчитать характеристики выбросов рассматриваемых процессов.

Исследуем периодически коррелированный случайный процесс [4] дискретного аргумента на ограниченном интервале в равноотстоящие моменты времени 11, ...,£^, Дt = ^+1 — с период ом Т = рД£, где р — такое, что Т < tN. Его математическое ожидание Е£(^), дисперсия (^) и корреляционная функция Щ(^,tj) удовлетворяют

соотношениям [5]

Е£(Л + Т) = Е£(ti), (ti + Т) = (Л), Щ(Л + T, tj + Т) = Щ(Л, tj)• (2)

На рис. 1 приведены средние М(г), дисперсии Б (г) и коэффициенты корреляции Я(г,]) для различных временных сдвигов, вычисленные по реальным данным (станция “Свердловск”) для января, г, ] = 1,..., 10. Видно, что характеристики реальных

-19 І

-15

-17

Температура

ІЛЛЛЛЛЛЛЛЛА

12345678910 Номер суток

Температура

7

1 2345678910 Номер суток

1,2

1

0.8

0.6

Коэф фициент корр еляции Г (І, })

— г а і+і)

-----г (і, І+2)

-----Г&І+3)

...... Г&І+4)

12345678910 Номер суток

Рис. 1. Математические ожидания (а), дисперсии (б) и коэффициенты корреляции (в) для различных временных сдвигов

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

1. Стохастическая модель периодически коррелированных временных рядов температуры воздуха

В основу численной вероятностной модели периодически коррелированных случайных процессов с периодом р положен алгоритм моделирования последовательности стационарно связанных гауссовских р-мериых векторов | = (£н,£2г, . . . ,Срг)Т, * = 1, . . . ,п, с ковариационной матрицей блочно-теплицева вида,

До Д1 ... Дп-1

Д( п) = ДТ До ... Дп-2

^ е 2 ... До

где Дк, к = 0,..., п — 1 — матрицы р х р. Вещественная положительно определенная матрица Д(п) симметрична, но ее блочные элементы могут быть несимметричными. Без нарушения общности мы можем рассматривать нормированную ковариационную матрицу Д(п),

Алгоритм моделирования последовательности векторов £*, * = 1,..., п,с корреляционной матрицей (3) основан на методе условных математических ожиданий [6], который реализуется по схеме

Iк В [к 1]<^(к—1)£(к-1) + Ск-1^к, к ^ . . . , п (4)

где ^1, ^2,..., — независимые гауссовские векторы размерности р, такие, что

Е^к^Т = 1р, Е^к^Т = 0, к = I. Здесь 0 — нулевая, а — единичная матрицы размерности р х р; 7(к) — перестановочная блочная матрица с к х к блоками размера р х р, на побочной блочной диагонали которой стоят блоки /р, остальные блоки — нулевые матрицы,

1(к) = (1Т,..., 1Т )т, в [к] = (вТ [к],...,ВТ [к])т, к = 1,...,п — 1,

В* [к] — матрицы раз мера р х р, а С* — нижние треугольные матрицы р х р, такие, что (к = СкСТ, где ( = До — ВТ [к] Д(к)В [к] — условная ковариационная матрица вектора |к ПРИ фиксированных значениях |к-1,..., |1,

Численная реализация схемы (4) требует вычисления на каждом шаге матрицы В [к] и остаточной ковариационной матрицы (к, В данной работе для этой цели используется рекуррентный алгоритм

ВТ [1] = ДТ Д-1, ВТ [1]= Д1Д-1, до = До, (о = До,

(ВТ[к +1],..., вТ[к +1]) = вТ[к] — вТ+1[к + 1]ВТ [к] /(к),

(ВТ[к + 1],... , Д[к + 1]) = ВТ[к] — ВТ+1[к + 1]ВТ [к] /(к), ( ,

Вк+1 [к + 1] = дк 1(Дк+1 — К*^(к)В [k]),

^к+1 [к +1] = д-1(ДТ+1 — иТ /(к)В [к]),

( = До — ИТ В [к], (5 к = До — ИТ В [к], Ск СТ = (, где Ик = (ДТ,..., ДТ )Т, И к = (Д1,..., Дк)Т, к = 1,... ,п — 1.

Составим из векторов &*, г = 1,..., п, построенных с помощью алгоритма (4), (5), вектор

&(п) = (£1Ъ £2Ъ . . . , СрЪ ^, . . . , £р2> . . . , £1п £2п . . . ,£рп)Т. (6)

Легко показать, что последовательность компонентов данного вектора является периодически коррелированной. Первые два условия из (2) обеспечиваются умножением каждого элемента последовательности (6) на элементы неслучайной периодической последовательности соответствующих стандартных отклонений с периодом Т = р и прибавлением к последовательности (6) неслучайной периодической последовательности соответствующих средних с тем же периодом. Отметим, что специфика периодической

Дк

Д(п)

Для моделирования бесконечных последовательностей

^1, £2, . . . , £р, £р+1, £р+2, . . . , £2р, . . . , £(п-1)р+1, £(п-1)р+2, . . . , £пр, . . .

может быть использована многомерная модель авторегрессии порядка п — 1 в виде & = ВТ [п — 1]&*-1 + ... + ВТ-1[п — 1]&-(п-1) + Сп-1^

где в качестве начальных векторов используются векторы &1,..., £п-1. Если выполняются условия стационарности процесса авторегрессии, то для него также выполняются соотношения (2),

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

2. Численная оценка характеристик выбросов температуры воздуха

Пусть событие В*(Ь, С) означает, что температура воздуха в г-й срок наблюдений опустилась ниже уровня С и на протяжении Ь суток в этот г-й час суток остается ниже данного уровня:

В*(Ь, С) = (0, < С, 0+м < С,..., < С}, г = 1,... ,р. (7)

Таким образом, мы исключаем разрывы серий, обусловленные существенными изменениями температуры в течение суток. Вероятности этих событий обозначим через р* = Р(В*(Ь, С)), г = 1,... ,р.

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

В качестве примера приведем зависимость от порядка модели авторегрессии вероятностей выхода температуры воздуха за уровень С = —20 °С в сторону более низких температур па протяжении интервалов времени длительностью от 20 до 30 дней (рис, 2), От порядка модели авторегрессии существенно зависит поведение ковариаций между значениями элементов ряда, удаленных друг от друга па большие временные расстояния, поэтому важны вопросы о достоверности “хвоста" ковариационной функции, оцененной по исходным данным, и вопрос о выборе порядка модели авторегрессии. Дня решения этой задачи представим модельные данные по аналогии с реальными. По каждой выборке £;Пт, 3 = 1,..., N гДе N = 104, сформированной с помощью модели, оценивается матричная ковариационная функция Щ0, Щ,..., Щп-1 (первая блочная строка матрицы Щп)) и по выборке автоковариационных функций при 3 = 1,..., 104

тк(т) = Шк,к],Ё{[к,к],...,Ё3т [k,k], к =1,...,p,

где Щ [к, к] — диагональные элементы блока Щ, вычисляются различные статистические характеристики этой оценки,

В таблице приводятся некоторые элементы модельной и фактической матричных корреляционных функций т\(т) и т1(т), полученные в результате нормировки соответствующих ковариационных функций, границы доверительного интервала тк(т) — ак, тк(т) + ак. Здесь ак — стандартное отклонение оценки корреляционной функции по выборке £п,т, оде т = 47 лет, В последней колонке таблицы приведено относительное стандартное отклонение ак/тк(т) оценки ковариационной функции, вычисленной по выборке ^Ппт, 3 = 1,..., 104, Из данных этой колонки видно, что относительная погрешность при временных сдвигах более девяти суток превышает 30%, а при сдвигах более двенадцати суток — 40%, При увеличении временного сдвига эта погрешность растет и может превышать величину 300%, Отсюда можно сделать вывод, что если

Рис. 2. Вероятности выхода температуры воздуха за уровень С = —20 °С в г-й срок длительностью Ь = 20, 21,... , 30 дней для различных порядков модели авторегрессии: а — п = 5, б — п = 15

Элементы матричной корреляционной функции, полученные в результате нормировки ковариантных функций

т Гк(т) - ак Гк{т) г*к(т) Гк(т) + а к Относительная ошибка СГк/Гк(т), %

0 1.0 1.0 1.0 1.0 0

1 0.736072 0.754693 0.755498 0.773314 2.46736090

2 0.501530 0.535056 0.536594 0.568582 6.26588618

3 0.360131 0.401991 0.404583 0.443851 10.41316845

4 0.278350 0.325671 0.328854 0.372992 14.53030819

5 0.175087 0.227822 0.231459 0.280558 23.14789616

6 0.126184 0.181028 0.185158 0.235872 30.29586583

7 0.121967 0.176435 0.180674 0.230904 30.87199252

9 0.126654 0.180399 0.184881 0.234144 29.79229375

10 0.108004 0.163093 0.167646 0.218183 33.77827375

12 0.081349 0.139407 0.148139 0.197465 41.64640226

15 0.012306 0.076859 0.144555 0.141412 83.98886272

20 -0.05461 0.019574 0.159629 0.093765 379.0283029

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

Таким образом, если мы хотим с заданным уровнем погрешности (менее 30 %) учесть начальные корреляции в моделируемом процессе, то порядок модели авторегрессии целесообразно выбирать не превосходящим 9,

В качестве второго способа выбора порядка модели рассмотрен подход, когда порядок модели выбирается из условия, что модельные вероятности выхода температуры воздуха за заданные уровни в течение различных заданных сроков близки к реальным. Все результаты получены для стационарного процесса и для первого срока измерения г = 1, т. е, для 0 часов. Рассматривались уровни С = Е£ь Е£х — 3, Е£х — 6, Е£х — 9 °С, где = — 15.2421 °С.

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

Для проверки адекватности модели реальным данным рассчитаны вероятности выхода температуры воздуха за заданные уровни С = 0, —5, —10, —15, —20, —25 °С длительностью Ь = 2 и Ь = 5 (дней) в зависимости от срока наблюдений для порядка модели авторегрессии п = 9, Число модельных реализаций N=10 000. Расчеты показали приемлемую точность модели (рис, 3),

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

СЬ

данного уровня:

В(Ь,С) = К; < С,6+1 < С, . . .,&+£„ < С}, (8)

Рис. 3. Вероятности выхода температуры воздуха за уровень С в 1-й срок наблюдения длительностью Ь = 5 (дней) по реальным (.......) и модельным (------) данным, п = 9

Температура

------ Уровень для модели ПКСП

------Уровень для 1 срока

-----Уровень для 2 срока

------ Уровень для 3 срока

------ Уровень для 4 срока

Номер суток

С

нарных процессов для различных сроков наблюдений, Р = 0.7, число модельных реализаций N=10000

где ] — номер начального элемента в событии Б(Ь,С), Требуется определить уровень С при заданной вероятности Р = Р(Б(Ь, С)) наступления события Б(Ь, С) и заданной длительности выхода Ь,

С

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

Список литературы

[1] Каган Р.Л., Капашкип В.К., Федорчепко Е.И. О расчете характеристик временных рядов методом статистическоі'о моделирования /7 Тр. Гл. геофиз. обсерватории им. А.И. Воейкова. 1972. Вып. 286. С. 71 82.

[21 Ф ЕДОРЧЕПКО Е.И. О суточном ходе характеристик выбросов температурных рядов /7 Тр. Гл. і'еофиз. обсерватории им. А.И. Воейкова. 1977. Вып. 397. С. 27 34.

[3] Гапдип Л.С., Каган Р.Л. Статистическая интерпретация метеорологических данных. Л.: Гидрометеоиздат, 1987. 320 с.

[4] Драган Я.П., Рожков В.А., Яворский И.Н. Методы вероятностного анализа ритмики океанологических процессов. Л.: Гидрометеоиздат, 1987. 320 с.

[5] Боков В.Н., Лопатухин Л.И., Микулинская С.М. и др. О межгодовой изменчивости волнения // Пробл. исследования и мат. моделирования ветрового волнения. СПб.: Гидрометеоиздат, 1995. С. 446-454.

[6] Ogorodnikov V. A., Prigarin S. М. Numerical Modeling of Random Processes and Fields: Algorithms and Applications, Utrecht: VSP. The Netherlands, 1996.

Поступила в редакцию 20 февраля 2008 г.