стохастическое моделирование и прогноз загрязнения атмосферы сиспользованием нелинейной регрессии
А. Пашкова
Содержание
ВВЕДЕНИЕ
1. ОБЗОР И ПОСТАНОВКА ЗАДАЧИ
2. МЕТОДОЛОГИЯ РЕШЕНИЯ
2.1 Модель
2.2 Алгоритм
3.3 Описание методов
2.4 Программа
3. ОПИСАНИЕ РЕЗУЛЬТАТОВ
3.1 Исходные данные
3.2 Подготовительный этап
3.3 Преобразование предикторов
3.4 Построение регрессионного уравнения
3.5 Оценка эффективностипрогностической схемы
3.6 Сравнительный анализ результатов
3.7 Результаты для других городов
ЗАКЛЮЧЕНИЕ
ЛИТЕРАТУРАВВЕДЕНИЕ
Разработка методов прогнозированиязагрязнения воздуха является одной из важнейших задач, возникающих в рамкахпроблематики охраны воздушного бассейна. Прогнозы и предупреждения о высокомуровне загрязнения воздуха служат основанием для проведения мероприятий порегулированию выбросов и уменьшению антропогенной нагрузки на окружающую средув периоды неблагоприятных метеорологических условий.
Актуальность указаннойтематики обусловлена тем, что в последние годы, несмотря на применяемыеприродоохранные меры, проблема чистоты атмосферы городов Российской Федерациине только не решена, но даже обострилась. Как следует из анализа результатовнаблюдений, за последнее десятилетие в крупнейших (численностью более 500 тыс.жителей) городах России, высокий уровень загрязнения воздуха сохранился и,согласно прогнозу, такая тенденция будет иметь место в течение ряда лет.Сравнение средних за год концентраций примесей с национальными стандартамикачества воздуха показало, что они превышают установленные нормативы предельнодопустимых концентраций (ПДК). Максимальные концентрации, превышающие ПДК вдесятки раз, регулярно регистрировались в большинстве (55-80%) крупнейшихгородов страны.
Проводимая в нашей странеприродоохранная политика предусматривает необходимость регулирования (т.е.сокращения) выбросов в периоды неблагоприятных метеорологических условий (НМУ).Чтобы обеспечить практическую возможность такого сокращения, соответствующиеорганы управления, предприятия и др. должны быть оповещены заранее онеобходимости своевременного сокращения выбросов вредных примесей.
В этой связи практическийинтерес представляют краткосрочные прогнозы — большей частью сзаблаговременностью от 1 до 3 суток. Их внедрение при обеспечении регулированиявыбросов позволяет в ближайшее время без значительных затрат улучшить состояниевоздушного бассейна и предотвратить появление высоких уровней загрязнениявоздуха.
Целью дипломной работыявляется разработка стохастической модели, которая позволяет производитькраткосрочный прогноз абсолютных уровней загрязнения воздуха на территориигородов методами регрессионного анализа. В соответствии с указанной целью вработе были поставлены следующие задачи:
- выбор метода и разработкастатистической модели прогноза загрязнения;
- выборвычислительных алгоритмов;
- программная реализацияматематической модели;
- построениесоответствующие этой модели прогностических схем загрязнения воздуха по даннымизмерений и оценка их эффективности на основе данных, полученных в несколькихгородах;
- проведение анализарезультатов, полученных на основе прогностических схем.
Оперативные методыпрогноза загрязнения воздуха изложены в действующем «Руководстве по прогнозу загрязнениявоздуха» РД.52.04.306-92 [5] и внедрены во всех управлениях Росгидромета. Работыпо прогнозированию загрязнения воздуха проводятся в 250 городах РоссийскойФедерации, предупреждения его возможного роста передаются более чем на 5000предприятий, на которых принимаются конкретные меры по снижению выбросов внеблагоприятные периоды.
Однако, методы прогноза,используемые в оперативной практике Росгидромета, не позволяют предсказыватьнаибольшие концентрации примесей, формирующиеся в воздухе отдельных районовпромышленных городов. А ведь предотвращение появления высоких концентрацийимеет существенное значение для решения проблемы защиты атмосферы отзагрязнения в период НМУ.
Модель, разработанная вданной работе, не имеет подобных недостатков, так как позволяет прогнозироватьне относительные характеристики загрязнения воздуха на территории городов, а ихабсолютные уровни.
1.ОБЗОР И ПОСТАНОВКА ЗАДАЧИ
Основная задача даннойработы ¾ краткосрочный прогноз максимальныхконцентраций загрязняющих примесей. Рабочая гипотеза состоит в том, чтоформирование опасного загрязнения воздуха на достаточно большой территориигорода обусловлено процессами рассеивания антропогенных выбросов, то есть, процессами,регулярно встречающихся при обычных технологических режимах работы промышленныхисточников или в суточном ходе изменения в интенсивности движенияавтотранспорта. Задача прогноза последствий аварийных выбросов вредных веществв атмосферу вследствие нарушения технологических процессов на предприятиях вданной работе не рассматривается.
Количественной оценкойнаибольшего загрязнения может служить максимальная (CМАХ) из измеренных за день концентраций на посту. Измеренные накаком-то конкретном посту концентрации значительно варьируются в течение сутокза счет изменчивости:
- выбросов,
- «внешних»метеорологических условий,
- локальныхфлуктуаций концентраций, связанных с флуктуациями локальных метеопараметров,определяющих перенос и рассеивание примесей.
При переходе к дневныммаксимумам влияние локальных флуктуаций сильно уменьшается, уменьшается такжероль систематических изменений выбросов в течение суток, так что связь с внешнимиметеорологическими условиями становится более устойчивой.
Решение основной задачи работысводится к разработке метода прогноза максимальной за день концентрации в точкенаблюдения. При такой постановки задачи необходимо выполнение следующеготребования: разрабатываемый метод должен быть эффективным как для области частонаблюдаемых значений, так и для высоких квантилей (экстремумов) функциираспределения CМ.
При наличии многолетнихрядов наблюдения за загрязнением атмосферы для этой цели можно применить хорошоразработанный математический аппарат регрессионного анализа. Он позволяет наоснове статистического анализа зависимостей между переменными прогнозироватьодну из них, если известны значения других.
Воспользуемся для решенияпоставленной задачи методом регрессионного анализа с определеннымипреобразованиями — достаточно простым и удобным математическим аппаратом.
В рамках регрессионногоанализа модель представляется в виде:
/>, (1)
где СМАХ –предиктант (в нашем случае максимальная концентрация рассматриваемой примеси засутки), xi – предикторы (в качестве предикторов используют различныеметеорологические характеристики), а bi – коэффициенты регрессии,которые требуется оценить.
e ¾ вектор ошибки (остатка). Предполагается, что e ¾ независимая случайная величина, имеющая нормальноераспределение N(0, s2).
Уравнение решаетсяметодом наименьших квадратов, т.е. из условия минимума среднего квадратаошибок. После того как определены параметры bi, получаем прогностическое уравнение,которое отличается от исходного модельного уравнения тем, что не содержатслучайной ошибки.
Успешное применениевыбранного математического аппарата требует обеспечить выполнение следующихусловий: двумерное (совместное) распределение плотности вероятности переменных(предиктанта с каждым из предикторов) подчиняется нормальному закону; формасвязи между переменными должна быть близкой к линейной. При этом, однако,заранее известно, что экстремумы случайной величины обычно не распределенынормально, и что связи между предикторами и предиктантом могут оказатьсянелинейными. Поэтому перед тем, как строить уравнение регрессии, необходимопреобразовать переменные таким образом, чтобы линеаризовать и нормализоватьсоответствующие связи.
После соответствующихпреобразований переменных, включающих нормализацию предиктанта ипредварительное исключение нелинейности связей, получаем рабочее уравнение:
/>, (2)
Задача с преобразованнымипредикторами решается методом пошагового регрессионного анализа. Данный виданализа позволяет включать в схему только те факторы, которые имеют значимуюкорреляцию с показателями загрязнения. Применение такого аппарата обусловленотем, что нет никакой гарантии, что между используемыми предикторами отсутствуеттесная корреляционная связь. Если же такая связь существует, то соответствующаясистема уравнений метода наименьших квадратов, используемая для определениякоэффициентов в уравнении регрессии, оказывается плохо обусловленной, а еерешение может привести к накоплению вычислительных ошибок.
После того как определеныпараметры bi, получаем прогностическое уравнение. По этому уравнениюрассчитываются прогностические значения максимальной концентрации загрязняющейпримеси. Применимость этого уравнения проверяется его испытанием на независимойвыборке.
Из значений СМАХи СМАХПРОГ, полученных с использованием прогностических уравнений позависимому и независимому рядам, формируется таблица результатов прогноза ирассчитываются статистические характеристики эффективности прогнозамаксимальной концентрации примеси.
Исходные данные дляразработки стохастических моделей были предоставлены ГУ «ГГО» по таким городам,как Санкт-Петербург, Обнинск, Милан, Мадрид, Новосибирск и др.
2.МЕТОДОЛОГИЯ РЕШЕНИЯ 2.1 Модель
Регрессионный анализ —это эффективный метод, который позволяет анализировать значительные объемыинформации с целью исследования вероятной взаимосвязи двух или большепеременных.
В регрессионном анализерассматривается связь между одной, зависимой, переменной и несколькими другиминезависимыми переменными. Эта связь представляется с помощью математической модели,то есть уравнением, которое связывает зависимую переменную с независимыми. Врамках регрессионного анализа модель представляется в виде:
/>, (3)
где СМАХ –предиктант (в нашем случае максимальная концентрация рассматриваемой примеси засутки), Xi – предикторы (в качестве предикторов используютразличные метеорологические характеристики и концентрации других загрязняющихпримесей), а bi – коэффициенты регрессии, которые требуется оценить.
Регрессионный анализиспользуется по двум причинам.
1. Описание зависимости междупредикторами и предиктантом помогает установить наличие возможной причиннойсвязи.
2. Получение аналитической зависимостимежду переменными дает возможность предсказывать будущие значения СМАХпо значениям предикторов.
Успешное применение этогоматематического аппарата требует выполнение двух условий:
1. Функции распределения переменных(предиктанта и каждого из предикторов) подчиняются нормальному случайномузакону.
2. Форма связи между переменными должнабыть близкой к линейной.2.2 Алгоритм
Предварительный этапразработки прогностической схемы состоит в подготовке исходного ряда данных:
1. Ряд разбивается на «обучающую» и«независимую» выборки. В данной работе прогностическая модель загрязнения атмосферыразрабатывается с использованием длительного ряда данных наблюдений. Две третиряда рассматриваются, как «обучающая» выборка для построения прогностической схемы,а оставшаяся одна треть применяется для проверки её эффективности нанезависимом материале (т.е. как «независимая» выборка). К «независимой выборке»относятся данные наблюдений, соответствующие неделям года с номерами, кратнымитрём (т. е. третья, шестая, девятая и т.д. недели). Остальные данные относятсяк «обучающей» выборке.
2. По «обучающей» выборке строитсяфункция распределения суточных максимумов концентраций и определяется её 60-ыйпроцентиль С60.
3. Устанавливается граничное значение СГРдля прогноза суточных максимумов, которое принимается равным С60.
Прогноз CMAX осуществляется по следующимправилам:
1. Если максимальная за предыдущие суткиC’MAX концентрация была ниже СГР,то прогнозируемая максимальная концентрация на очередные сутки CMAXПРОГ принимается равной C’MAX(«инерционный прогноз»).
2. Если максимальная за предыдущие суткиC’MAX концентрация была выше или равна СГР,то прогноз осуществляется с использованием прогностических схем.
Применение методалинейной регрессии требует, чтобы корреляционные связи между предиктантом скаждым из предикторов были близки к линейным, однако это условие не всегдавыполняется. Для исключения нелинейности связей предикторы нужно преобразоватьс помощью кривых зависимости показателя загрязнения воздуха от отдельныхметеопараметров, построенных по использованному для разработок материалунаблюдений. При этом каждое значение предиктора меняется на соответствующее емусреднее значение характеристики загрязнения.
- Для каждойградации предиктора (их должно быть не менее 5) рассчитать среднее значение />CMAX. При недостаточном количествеслучаев в одной из градаций, она объединяется с соседних. Таким образом,получаем набор точек с абсциссами M(CMAX) и ординатами, соответствующими серединам отрезковосреднения.
- Построить графиккусочно-линейной функции, у которой полученные точки являются угловыми.
- Каждому значениюпреобразованного предиктора сопоставляется значение кусочно-линейной функции всоответствующей точке.
Связь преобразованныхтаким образом предикторов с предиктантом в значительной степени линеаризуется.Этот прием позволяет учесть реальный вид связи в каждом конкретном случае. Онблизок к так называемому «кусочно-линейному» преобразованию, применяемому припостроении моделей для прогноза погоды.
При возникновениитрудностей, связанных с тем, что данные, подчиняющиеся какому-нибудьнесимметричному распределению, должны быть подвергнуты анализу, теория которогоразработана в основном для нормального распределения, можно преобразоватьэмпирическое распределение в нормальное («нормализовать переменные») и затемпродолжить анализ на базе известной теории.
Для нормализациипеременных используется стандартное преобразование выборочной функциираспределения в нормальную (гауссову) со средним, равным 0, и стандартнымотклонением, равным 1. Это преобразование осуществляется по формуле
/>, (4)
где Ф-1(t) –обратная функция к функции распределения нормальной случайной величины сосредним значением ноль и стандартным отклонением единица, а F(x) — выборочнаяфункция распределения рассматриваемой случайной величины X.
Задача с преобразованнымипредикторами решается методом многомерной пошаговой регрессии. На каждойитерации этого метода ищется предиктор, имеющий наибольшую связь спредиктантом. Таким образом определяются наиболее значимые предикторы, которыеследует включить в уравнение регрессии.
Если значимыми оказалисьдва предиктора, соответствующие двум срокам измерения одного и того жеметеорологического параметра, то в уравнение регрессии включается тот, которыйбольше связан с предиктантом. В итоге должны остаться 4 – 7 наиболееинформативных предикторов, связь которых с предиктантом наиболее значима.
Данный вид анализапозволяет включать в схему только те факторы, которые имеют значимую корреляциюс показателями загрязнения. Применение такого аппарата также обусловлено тем,что нет никакой гарантии, что между используемыми предикторами отсутствуеттесная корреляционная связь. Если же такая связь существует, то соответствующаясистема уравнений метода наименьших квадратов, используемая для определениякоэффициентов в уравнении регрессии, оказывается плохо обусловленной, а ее решениеможет привести к накоплению вычислительных ошибок. После того как определеныпараметры bi, получаемстохастическую модель процесса, которая может быть для краткости представлена ввиде:
/>, (5)
Здесь [Xi] –преобразованныепредикторы, I – количество использованныхпредикторов, b0и bi – коэффициенты регрессии. Значения b0и bi определяетсяс помощью метода наименьших квадратов. По этому уравнению рассчитываютсяпрогностические значения максимальной концентрации загрязняющей примеси.
Из значений Смах и Смахпрог, полученных сиспользованием прогностических уравнений по зависимому и независимому рядам,формируется таблица результатов прогноза. Рассчитываются статистическиехарактеристики эффективности прогноза максимальной концентрации примеси.
Эффективностьразработанных прогностических схем проверяется по зависимым (использованным дляпостроения уравнения регрессии) и независимым (не использовавшимся для построенияуравнения регрессии) данным наблюдений./>
Оправдываемостьиндивидуального прогноза максимальной концентрации примеси Смахпрог за конкретные суткиоценивается путемсопоставления этой прогностической концентрации с определенной по даннымнаблюдений фактической максимальной за сутки концентрацией Смах. Прогноз считаетсяоправдавшимся, если при Смах> ПДК выполняется условие:
/>, (6)
или если при Смах/> ПДК выполняется условие
/>, (7)
где ПДК — установленнаяМинздравом РФ максимальная разовая предельно допустимая концентрация примеси ватмосферном воздухе населенных мест. 2.3 Описание методов
1. Нахождение СГР.
n-й процентиль — этотакое значение, ниже которого расположено n процентов наблюдений рассматриваемойпеременной. График функции распределения случайной величины X имеет ступенчатый вид. Значениефункции F(X) равно:
/>, k = 0…M-1, (8)
где M – объём выборки, а k – порядковый номер события вупорядоченном по возрастанию массиве. Как известно, то α-квантильоднозначно задаётся уравнением: F(xα) =α. Значит за 60 процентиль можно принять элемент с порядковым номером k = 0.6M (округление производим в большую сторону).
2. Нормализации.
Нормализацияосуществляется по формуле:
/>, (9)
График функциираспределения случайной величины Xимеет ступенчатый вид. Значение функции F(X) равно:
/>, k = 0…M-1, (10)
Так как при k = 0 F(Xk) обращается в ноль, то [Xk] становится равным минус бесконечности, что является нежелательным,заменим формулу (10) на:
/>, k = 0…M-1.(11)
При достаточно больших Mпогрешность в значениях F(Xk), вычисляемых по формуле(11) становитсямала. При этом F(X) нигде не обращается в ноль или M, а значит [Xk] принимают только конечные значения.
Вместо функции, обратнойк функции распределения нормальной случайной величины, Ф-1 можноиспользовать её аппроксимацию (погрешность e-16).
3. Пошаговая регрессия.
Имеется набор независимыхпеременных X1…Xn, которые являются кандидатами на роль предикторов СМАХ, ислучайная выборка объема М. Рассмотрим стандартную пошаговую процедуру(F-метод), которая состоит из правила включения переменных и правилаисключения. Включение и удаление переменных осуществляются с помощью критерия,который имеет F-распределение, и называется либо F-включения, либо F-удаления.
Более точно, предположим,что в набор с уже включено k переменных, k = 0, 1… M-1. Тогда значение F-включения дляпеременной X (не входящей в с) вычисляется по формуле:
/>, (12)
где rСмахX*с – множественный коэффициент корреляции. Величина Fслужит статистикой критерия для проверки гипотезы о том, что предсказание СМАХзначимо не улучшается при включении X в набор с. Аналогично, величинаF-удаления для какой-либо переменной X из с служит статистикой критериядля проверки гипотезы о том, что набор с' получающийся из с приудалении X и содержащий k'=k-1 переменных, предсказывает СМАХ так жехорошо как и набор с.
/>, (13)
Правило остановки, обычноиспользуемое в стандартной процедуре, основано на задании допустимого минимумаF-включения. По умолчанию предполагается, что минимум F-включения равен 4. Дляудаляемых переменных также выбирается допустимый минимум. F-удаления величинадолжна быть меньше минимума F-включения. По умолчанию принимается, что минимумF-удаления равен 3,9. Рассмотрим теперь подробно шаги стандартной процедуры.
Шаг 0. Вычисляются величины F-включениядля i = 1…n. Статистика критерия даетсявыражением
/>, (14)
которое получается изформулы (?) подстановкой k=0.
Шаг 1. Переменная Xi1, которойотвечает наибольшее значение F-включения (или, что эквивалентно, наибольшаявеличина квадрата коэффициента корреляции с СМАХ), выбирается какнаилучший предиктор. Величина F-удаления для Xi1 в этом случаесовпадает с величиной F-включения. Далее вычисляются значение F-включения иF-удаления.
Если все вычисленныезначения F-включения меньше установленного минимума, то далее выполняется шагS.В противном случае происходит переход на шаг 2.
Шаг 2. Переменная Хi2, имеющаянаибольшее значение F-включения (или, что эквивалентно, наибольший квадрат частногокоэффициента корреляции с СМАХ при фиксированном значении Xi1),выбирается как наилучший предиктор для СМАХ при условии, что ужевыбрана переменная Xi1. Если все значения F-включения меньшеустановленного минимума, то далее выполняется шаг S. В противном случаепроисходит переход на шаг 3.
Шаг 3. а) Пусть L обозначает набориз l независимых переменных, которые включены в уравнение регрессии.Если какое-либо из значений F-удаления для переменных из L меньше, чем соответствующийминимум, то переменная, которой соответствует наименьшее значение F-удаления,удаляется из набора и выполняется шаг Зb с заменой l на l-1. Еслидля всех переменных, не входящих в L, значение F-включения меньшеустановленного минимума, то выполняется шаг S. В противном случае в набор Lдобавляется переменная, которой соответствует максимальное значениеF-включения, и l заменяется на l + 1. b) Вычисляются значения F-удаления между СМАХи переменной Хi из L при заданных остальных l-1 переменных изL и значение F-включения между СМАХ и каждой переменной Xi,не входящей в L, при данных переменных из L.
Шаги 4, 5... Рекуррентно повторяется шаг 3. шаг Sвыполняется а) если F-включения для всех переменных, не входящих в L, меньшеустановленного минимума, b) если для всех переменных из L значение F-удалениябольше установленного минимума или с) число включенных переменных равно р.
Шаг S. Вывод результатов.
4. Метод наименьшихквадратов. Методнаименьших квадратов — простой и быстрый способ получить неизвестные параметрыв функциональной зависимости (/>) и оценить их погрешности. Минимизируетсясумма квадратов отклонений реально наблюдаемых CMAX от их оценок CMAXПРОГ (имеются в виду оценки с помощьюпрямой линии, претендующей на то, чтобы представлять искомую регрессионнуюзависимость):
/>, (15)
Для решения задачирегрессионного анализа методом наименьших квадратов вводится понятие функцииневязки:
/>, (16)
Условие минимума функцииневязки – равенство нулю всех частых производных, то есть:
/>, (17)
Полученная системаявляется системой I+1 линейныхуравнений с I+1 неизвестными b0...bI. Так как матрица системы уравненийявляется симметричной, то система решается методом LU-разложения.
5. Сортировка. Сортировка массивов осуществляетсяметодом />разделения (Quicksort). Шаги алгоритма таковы:
1. Выбираем вмассиве некоторый элемент, который будем называть опорным элементом. С точкизрения корректности алгоритма выбор опорного элемента безразличен. С точкизрения повышения эффективности алгоритма выбираться должна медиана, но бездополнительных сведений о сортируемых данных её обычно невозможно получить.
2. Операцияразделения массива: реорганизуем массив таким образом, чтобы все элементы,меньшие или равные опорному элементу, оказались слева от него, а все элементы,большие опорного — справа от него. Обычный алгоритм операции:
— два индекса — lи r, приравниваются к минимальному имаксимальному индексу разделяемого массива соответственно;
— вычисляетсяопорный элемент m;
— индекс lпоследовательно увеличивается до m или до тех пор, пока l-й элемент непревысит опорный;
— индекс rпоследовательно уменьшается до m или до тех пор, пока r-й элемент не окажется меньше опорного;
— если r = l — найдена серединамассива — операция разделения закончена, оба индекса указывают на опорныйэлемент;
— если ll и r, которые были достигнуты. Следуетучесть, что если какая-либо граница (l или r) дошла до опорного элемента, то при обмене значение mизменяется на r или l соответственно.
3. Рекурсивно упорядочиваем подмассивы,лежащие слева и справа от опорного элемента.
4. Базой рекурсииявляются наборы, состоящие из одного или двух элементов. Первый возвращается висходном виде, во втором, при необходимости, сортировка сводится к перестановкедвух элементов. Все такие отрезки уже упорядочены в процессе разделения.
Алгоритм называетсябыстрой сортировкой, поскольку для него оценкой числа сравнений и обменовявляется O(n*lg n).
2.4 Программа
Для автоматизации расчётапо описанной модели была составлена прикладная программа. Она была написана наязыке C++ в среде разработки Microsoft Visual Studio 2005, и её объёмсоставляет несколько тысяч операторов. Программа состоит из трёх основныхмодулей: управляющая часть, блок с функциями, реализующими логику всехнеобходимых методов, и модулем, который отвечает за взаимодействие программы спользователем. Таким образом, программа составлена в соответствии с шаблономпроектирования модель-вид-контроллер (model-view-controller).
Управляющая часть –контроллер – позволяет разграничить модель и представление модели. Она включаетв себя последовательность команд, которые вызывают требуемые функции.
Основные функцииреализованы во втором модуле программы. В нём реализованы все этапы,необходимые для построения прогностической схемы, начиная от разбиения всехданных на зависимую и независимую выборки, и заканчивая прогнозированиемконцентрации загрязняющего вещества для введённых параметров.
Основные функции,реализованные во втором модуле программы:
- Считывание данныхиз файла. Для удобства пользователя исходные данные читаются из файла в форматеxls.
- Деление выборкина зависимую и независимую части. В «обучающую» выборку попадают первые двенедели из каждых трёх, а в независимую – каждая третья неделя. Отсчётначинается от первого понедельника года.
- Определениеграничного значения и отбрасывание из выборки неподходящих данных. В качествеаргумента данной функции передаётся процентиль, по которому определяетсяграничное значение. Если значение концентрации меньше граничного, тосоответствующая строка отбрасывается из выборки.
- Линеаризация.Функция преобразует предикторы так, чтобы исключить нелинейность связей. Дляэтого все значения предиктора делятся на несколько интервалов и для каждогоинтервала находится среднее значение концентрации. Таким образом, строитсякусочно-линейная функция, которая проходит через точки с координатами: сединаинтервала, соответствующее среднее значение. Координаты этих точек запоминаютсяв специальной структуре для дальнейшего использования.
- Нормализация.Данная функция для каждой преобразуемой величины вычисляет значение выборочнойфункции распределения вероятности. Затем она сопоставляет этому значениюзначение функции, обратной к функции распределения нормальной случайнойвеличины. Таким образом, происходит преобразование всех случайных величин кнормальным.
- Определениепредикторов, которые будут включены в модель. Необходимо выявить наиболеезначимые предикторы. Для этого используется метод пошагового регрессионногоанализа. При пошаговом алгоритме для критерия Фишера назначаются значениявключения и исключения переменных. За счёт этого можно регулировать количествовключаемых в модель предикторов.
- Функциясоставления системы для метода наименьших квадратов и её решения — определениекоэффициентов регрессии.
- Проверкаполученной схемы на «обучающей» и «независимой» выборках. Функция с помощьюсохранённых параметров для преобразования линеаризует, нормализует необходимыевеличины и вычисляет прогностическое значение концентрации загрязняющеговещества.
- Сохранение.Функция записывает все промежуточные данные и результаты вычислений в выходнойфайл out.xls. Это позволяет пользователю легко анализировать результаты,строить все необходимые графики и оценки для модели средствами MS Excel,OpenOffice Calc, Statistica и другие.
Также в этом модулепрограммы реализованы вспомогательные функции такие как:
- Вычислениефункции вероятности нормального распределения со средним 0 и ско 1.
- Вычислениеобратной функции вероятности нормального распределения.
- Вычислениеопределителя матрицы.
- Решение системы спомощью LU-разложения.
Благодаря третьемумодулю, который отвечает за визуализацию, пользователь имеет возможность,получать некоторые промежуточные результаты, в зависимости от них вводитьразличные параметры и корректировать работу программы.
Весь алгоритм программыможно представить в виде блок-схемы (рис.1).
/>
Рис. 1.Блок-схема алгоритма.
3.ОПИСАНИЕ РЕЗУЛЬТАТОВ 3.1 Исходные данные
Исходные данные (рис. 2)для разработки стохастической модели были предоставлены ГУ «ГГО» постанции, расположенной на ул. Пестеля (г. Санкт-Петербург). Эти данныехарактеризуют загрязнение атмосферного воздуха озоном за 2002 год.
В разрабатываемойстохастической модели связь межу предиктантом и предикторами описывается в виде/>,где
- /> — предиктант,максимальная за сутки концентрация озона (мкг/м3);
- /> – линейнаяфункция от n предикторов;
- в качествепредикторов /> используются следующие величины:
/> - максимальнаяконцентрация озона (мкг/м3) за предыдущие сутки;
/>-7- концентрация оксидаазота (мкг/м), измеренная в 7 часов;
/>-7 — концентрациядиоксида азота (мкг/м3), измеренная в 7 часов;
/> - концентрация озона(мкг/м3), измеренная в 7 часов;
/> - скорость ветра (м/с)в 6 и 12 часов;
/> - направление ветра(дес.град) в 6 и 12 часов;
/> - атмосферное давление(мб) в 6 и 12 часов;
/> - температура воздуха(°С) в 6 и 12 часов;
/> - относительнаявлажность воздуха (%) в 6 и 12 часов;
/> - атмосферные явления(шифр), наблюдаемые в 6 и 12 часов.
Длина массива данныхсоставляет 273.
/>
Рис.2.Исходные данные3.2 Подготовительныйэтап
Чтобы вычислитькоэффициенты функции /> при каждом из предикторовнеобходимо подготовить исходные данные. С этой целью:
1. Ряд значений /> максимальныхза сутки концентраций озона разбивается на обучающую и независимую выборки.Далее все преобразования производятся только с обучающей выборкой.
2. Осуществляетсяцензурирование выборки: сортируется массив данных в соответствие с ростомпеременной />, где /> - значение /> засутки до срока, на который дается прогноз. Находим значение 60%-ного квантиляфункции распределения />. — В рассматриваемом примерезначению 60%-ного квантиля функции распределения /> соответствуетконцентрация /> =/>63,42 мкг/м3. Делим исходные данные на две группы. Впервую группу, объем которой составляет 60% от общей выборки, войдут значения />/>, сопутствующие метеоусловия ирасчетные параметры, которые наблюдались при /> осуществляется по уравнению />. Вовторую группу, объем которой составляет 40% от общей выборки, — значения />, /> , />, />,сопутствующие метеоусловия и расчетные параметры, которые наблюдались при /> >=63,42 мкг/м3. Для этой группы прогноз /> осуществляется с использованиемпрогностических схем.
3.3 Преобразованиепредикторов
1. Линеаризация(рис.3). Исключение нелинейности связей между предиктантом /> ипредикторами.
/> />
/> />
Рис. 3. Линеаризация предикторов
2. Для повышенияэффективности модели в области высоких значений необходимо нормализовать всепеременные, входящие в модель (рис. 4). Для этого каждое значение переменнойсопоставляем значение обратной функции к функции распределения нормальнойслучайной величины со средним значением ноль и стандартным отклонением единица.Ее аргументом является значение выборочной функции распределения рассматриваемойслучайной величины.
/> />
Рис. 4. Графики функции распределенияпредиктанта (до и после преобразования)
3.4 Построениерегрессионного уравнения
Подготовленные исходныеданные применяются для получения коэффициентов уравнения. Применяем методпошаговой многомерной регрессии. Устанавливаем параметры: F-включения = 2,F-удаления = 1,9. После девяти итераций получаем предикторы, которые необходимовключить в систему. Так как два предиктора, F-6 и F-12, соответствуют одному итому же метеорологическому параметру, в систему включаем только тот, которыйбольше коррелирует с предиктантом, то есть F-6.
Таблица 1.Оценка предикторовF-включения F-удаления In/Out Предиктор 89,0213 In
/> 31,3170 In С7 3,4080 In F-6 6,1487 In D-6 4,7916 In NO2-7 2,2418 In F12 3,5270 In V12 3,0150 In Т-6 2,3492 In NO-7 0,4037 Out Р6 0,0001 Out Т-12 0,0001 Out D12 0,4391 Out V-6 0,0083 Out АЯ-12 0,9863 Out АЯ-6 0,6304 Out Р12
После того как определенынаиболее значимые предикторы, ищем коэффициенты b0и bi. Для этого используем метод наименьших квадратов. В результате получаем прогностическоеуравнение:
/>/>, (18)
Для того чтобы получитьокончательные значения концентрации />, необходимо выполнитьпреобразование, обратное нормализации. 3.5 Оценка эффективности прогностической схемы
По полученному уравнениюопределяется ожидаемое значение />.
Убедимся в эффективностипостроенной прогностической схемы с помощью корреляционного графика междуспрогнозированными и измеренными концентрациями озона для «обучающей выборки»(рис. 5).
/>
Рис.5.Корреляционный график измеренных и прогностических суточных максимумовконцентрации озона
Квадрат коэффициентакорреляции между данными измерения и прогноза составил 0,83, а угловойкоэффициент регрессии оказался равным 1,06.
Проверим применимостьпрогностической схемы на «независимой выборке». Корреляционный график,показывающий связь между прогностическими и фактическими концентрациями озона,приведен на рисунке (рис. 6).
/>
Рис.6.Корреляционный график измеренных и прогностических суточных максимумовконцентрации озона
Эффективность построеннойпрогностической схемы подтверждают выполненные оценки: квадрат коэффициентакорреляции между данными измерения и прогноза составил 0,78, а угловойкоэффициент регрессии оказался равным 1,09.
Оправдываемостьиндивидуального прогноза максимальной концентрации примеси Смахпрог законкретные сутки оценивается с помощью формул 6-7. Прогноз оказался оправдавшимсядля «обучающей выборки» в 89% случаев, а для «независимой выборки» — в 87%.
Таким образом, построенную схемупрогноза концентраций озона следует считать эффективной.
3.6 Сравнительныйанализ результатов
Сравним результаты, построенныес помощью приведённой выше прогностической схемы и с помощью методаинерционного прогноза. Полученные оценки сведём в таблицу. Прогностическая схема Инерционный прогноз «Обучающая выборка» «Независимая выборка» «Обучающая выборка» «Независимая выборка» Уравнение линии тренда
/>
/>
/>
/>
Квадрат коэффициента регрессии R2 0,8256 0,7808 0,7305 0,6585 Оправдываемость 89% 87% 75% 79% 3.7 Результатыдля других городов
/> />
Рис.7. Корреляционныйграфик измеренных и прогностических суточных максимумов концентрации озона(Милан)
прогнозконцентрация примесь программа
/>/>
Рис.8. Корреляционныйграфик измеренных и прогностических суточных максимумов концентрации озона (Winfield)
/> />
Рис.9. Корреляционныйграфик измеренных и прогностических суточных максимумов концентрации озона(Обнинск)
/> />
Рис.9. Корреляционный графикизмеренных и прогностических суточных максимумов концентрации озона(Новосибирск)
ЗАКЛЮЧЕНИЕ
При выполнении даннойработы получены следующие основные результаты:
1. Разработанстатистический метод прогноза максимальных за день концентраций примесей вотдельных точках города и построены соответствующие прогностические схемы (дляСанкт-Петербурга, Новосибирска, Милана). Проведенный анализ результатов играфические материалы показали, что применение выбранной стохастической моделипозволяет прогнозировать максимальные концентрации в период повышенногозагрязнения воздуха достаточно эффективно. Оценка эффективности схем поиспользованному и независимому материалу показала, что оправдываемость прогнозовнаибольших концентраций составила в среднем 80%, а их предсказуемость — 85%.
2. Разработана иотлажена компьютерная программа на алгоритмическом языке C++, реализующая указанный метод ипозволяющая прогнозировать суточные максимумы концентрации вредных примесей вгородах. С учетом универсальности использованного алгоритма, модно надеяться,что данная программа может быть также использована для решения задачипрогнозирования в других областях.
3. Показано, чтоприменение метода множественной линейной регрессии с предварительнымисключением нелинейности связей и нормализацией предиктанта позволяет успешнопрогнозировать максимальные концентрации и их экстремумы в период повышенногозагрязнения воздуха в городе. Необходимость преобразований переменных вызвананелинейной формой зависимости предиктанта от предикторов и асимметриейраспределения функции плотности вероятности.
Выполненная работанаправлена на повышение качества оценки уровня загрязнения на основеиспользования новой методики прогноза в целях повышения эффективности охранычистоты воздушного бассейна. Полученная модель может быть рекомендована дляоперативного использования в промышленных городах, в том числе для составленияпредупреждений об опасных уровнях загрязнения воздуха.
ЛИТЕРАТУРА
1. Афифи А., ЭйзенС. Статистический анализ: Подход с использованием ЭВМ. Пер. с английского. — М.: Мир, 1982. — 488с.
2. Берлянд М.Е.Современные проблемы атмосферной диффузии и загрязнения атмосферы.// Л.: Гидрометеоиздат, 1975.—448с.
3. Кириллова В.И. ….Автореферат диссертации …… кандидата геогр. наук
4. Норман Дрейпер,Гарри Смит. Прикладной регрессионный анализ. Множественная регрессия = AppliedRegression Analysis. — 3-е изд. — М.: «Диалектика», 2007. — С. 912.
5. «Руководство попрогнозу загрязнения воздуха» РД 52.04.306-92. // СПб.: Гидрометеоиздат, 1993.- 104 с.
6. Себер Дж.Линейный регрессионный анализ. М.: Мир, 1980, с. 456.