УДК 502.1
ЭКОЛОГИЧЕСКИЙ МОНИТОРИНГ И МЕТОДЫ МНОГОМЕРНОЙ МАТЕМАТИЧЕСКОЙ СТАТИСТИКИ Юрий Александрович Пичугин Российский государственный педагогический университет им. А.И. Герцена
yury-pi chu gi n@rambl er.ru
статистическое районирование, информативность, прямые и косвенные данные
Рассмотрены методы, облегчающие анализ многомерных массивов информации экологического мониторинга: статистическое районирование с использованием
выборочных главных компонент; отбор наиболее информативных данных прямых и косвенных измерений с использованием информации по Шеннону.
ENVIRONMENTAL QUALITY MONITORING AND MULTIDIMENSIONAL STATISTICAL ANALYSIS Yury A. Pichugin Herzen State Pedagogical University of Russia
yury-pi chu gi n@rambl er. ru
statistical division into districts, information, direct and indirect data
Methods for multidimensional environmental quality monitoring data analysis are considered: the statistical division into districts by using principal component analysis; the selection of direct and indirect measuring data by using Shannon information estimations.
Специфика экологического мониторинга многих районов нашей страны заключается в том, что процесс сбора информации нередко оказывается весьма непростым и трудоемким. В особенности это касается, например, районов Крайнего Севера. Указанное обстоятельство, в свою очередь, затрудняет синхронизацию во времени данных с различных точек (пунктов наблюдения), без которой невозможно делать объективные выводы об обстановке в том или ином регионе в целом и о динамике происходящих процессов. Наиболее эффективными методами преодоления указанной проблемы являются статистическое районирование территорий (выявление статистически однородных зон) и определение наиболее информативных точек сбора информации.
Статистическое районирование
Задача объективного районирования может быть решена методами многомерной математической статистики, в частности методом снижения размерности [2]. Таким образом, под объективным районированием мы будем понимать статистическое районирование.
Предположим, что по некоторому экологическому показателю имеются синхронизированные данные, которые можно представить в виде матрицы A={aj}, i=1,2,...m, j=1,2,...n (т.е. таблицы), где i - номер строки и пункта сбора информации, а j -номер столбца и порядковый номер по времени. Если в данных наблюдений имеют пропуски, то приходится удалять весь столбец, или заполнять пропуски, используя
101
пространственно-временную интерполяцию. От матрицы А необходимо перейти к
~ ^ ~ Он —
матрице А = {~} центрированных и нормированных значений Оц =-----------------------, где
1 " а/
1 П 9 1
2 Оц, а2 = — 2 (Он — а )2 — оценки среднего и дисперсии для каждого пункта. П]=1 1 П]=1 1
Далее следует вычислить матрицу взаимных корреляций К = — АА (' — знак
п
транспонирования) и ортогональную матрицу Р, приводящую Я к диагональному виду,
т.е. P'RP=K=diag(k{k2.Xm) так, чтобы значения X не возрастали с ростом номера i (diag
— стандартное обозначение диагональной матрицы). Вычисление матрицы Р осуществляется известными методами вращений с последующей перестановкой столбцов в порядке убывания значений X или итераций фон Мизеса, которые сразу обеспечивают
нужный порядок столбцов. Для проверки статистической гипотезы [1]
Н0:Xl>X2>.>Xh>Xh+l=Xh+2=.=Xm (т.е. для определения значения И) достаточно простым и надежным средством является правило "сломанной трости" [8,9]: в качестве И берётся
л 1/ + 1/ + + 1/
х к /к + 7 (1+к) + ••• + /ш максимальное значение, удовлетворяющее неравенству ------->---------------------------------------------------------, где
1гЯ ш
1тЯ — след матрицы Я, т.е. сумма элементов главной диагонали. Учитывая, что
ш ш
Ш = £ Гц =ХХ/ = ш, т.к. Я не ковариационная, а корреляционная матрица (данные
I=1 I=1
нормировались), следует проверять условие ХИ> ^^ + •• • + .
Обозначим через РИ матрицу, составленную из первых И столбцов Р. Каждая строка матрицы РИ соответствует конкретному пункту сбора информации. Стоки с одинаковой сигнатурой элементов соответствуют пунктам из одной однородной зоны, которые синхронно откликаются на изменение ситуации в исследуемом регионе в целом. В качестве примера использования вышеописанного метода можно предложить статью [7], где проведено статистическое районирование территории Танзании. В работе [7] столбцы матрицы РИ по традиции метеорологии и климатологии называются эмпирическими ортогональными функциями (ЭОФ), хотя общепринятое название в математической статистике - базис выборочных главных компонент [2].
Следует отметить, что в этой задаче совсем не обязательно использовать все столбцы матрицы РИ, что может привести к слишком большому количеству зон. При выборе этого числа первых 5 столбцов РИ (соответственно, будем иметь не РИ, а Рз) можно и следует руководствоваться соображениями иного характера (экономического, демографического и т.п.), соблюдая лишь условие з<И. Это ничуть не лишает указанный метод объективности, которая заключена в вычислении столбцов матрицы Р и значения И.
Оценка информативности пунктов наблюдения и показателей
Матрица РИ может быть использована и для решения другой задачи - задачи определения (выбора) наиболее информативных пунктов наблюдения. Столбец матрицы
А, содержащий полный набор сведений по некоторому показателю со всех пунктов наблюдения за один и тот же срок, обозначим через у, т.е. у имеет размерность т (по количеству пунктов сбора информации). Напомним, что А содержит центрированные
102
(отклонения от среднего значения) и нормированные данные. Гипотеза Н 0 (см. выше) или обоснованный соответствующим критерием выбор к позволяет предполагать, что в регрессионной модели вида [2,6] у=Ркг+є компоненты вектора погрешностей этой модели є взаимно распределены нормально, взаимно независимы и имеют одинаковую дисперсию. Компоненты вектора г, который имеет размерность к, и есть выборочные главные компоненты. Эти компоненты можно интерпретировать как скрытые факторы, определяющие обстановку по анализируемому показателю в целом, а соответствующие столбцы Рк или ЭОФ есть формы колебаний, подчиненных этим скрытым, неявным факторам.
Пусть 1,2, ^ подмножество пунктов наблюдения состоит из і элементов (1<і<да), тогда количество информации по Шеннону, содержащееся в компонентах вектора у с номерами из подмножества q относительно скрытых факторов (компонент вектора г), равно [5] Is(yq,z)=/2logЬdet(I+ЛhМ(q)). Здесь М^) - нормированная
информационная матрица, равная М(а)=о~2Р'к(а)Рк(о); Рн(о) содержит строки Рн с номерами из q; Лh=diag(k{k2...'khУ; выбор основания логарифма Ь (Ь>1) задает масштаб единиц информации и не имеет принципиального значения (не влияет на результат отбора); о2 - смещенная оценка дисперсии компонент вектора є вычисляется по формуле
оценивании [3]).
Начиная с =0 (т.е. q=0 - пустое множество) можно осуществить
последовательный отбор компонент у, увеличивая на каждом шаге отбора количество элементов q (7:=7+1) и максимизируя 1з(у^). Это приведет к упорядочиванию всех элементов по информативности, т.к. каждый следующий в порядке отбора элемент вносит не большее количество информации, чем предшествующий из числа отобранных ранее. Если все компоненты имеют нормальное распределение, то при общем росте количества информации в процессе последовательного отбора приращения количества информации монотонно убывают. В качестве ограничительного условия отбора наиболее информативных данных можно взять неравенство 1ц(у^)11з(у^)<Ь<1, где 5x100% — наперед заданный процентный уровень относительного количества информации; Is(y,z)=У2logbdet(I+ЛhМ) — полное количество информации относительно скрытых факторов; М=а~2Р'ьРк-
Поясним суть информативного отбора данных. Пусть отобрано по максимальной информативности I пунктов наблюдения (<т). Тогда вектору у нормированных и центрированных данных наблюдений, который имеет размерность т, соответствует вектор размерности I, обозначаемый yq и содержащий только те компоненты у, номера которых отобраны, т.е. входят в q. По данным, содержащимся в у^ можно восстановить остальные данные (компоненты у), пользуясь формулой уг=Рн(г)(Р'1,(о)Рк(у))1Р%Ы)у^ где Рй(г) содержит все строки Рй, номера которых не входят в q, т.е. г={1,2,..,,т}^ . Если на каждом шаге последовательного отбора выбор был безальтернативен (не было разных пунктов, дающих одинаковое приращение количества информации), то указанное выше восстановление недостающих данных имеет минимальную среднеквадратичную ошибку для фиксированного значения t.
Наглядным примером использования этого метода может служить работа [4], где определены наиболее информативные для метеопрогнозов по динамическим моделям зоны измерения начальных данных по Северному полушарию, т.е. зоны, где случайные погрешности метеорологических измерений могут приводить к существенным ошибкам прогнозов с использованием динамических моделей.
т
т - ^ і=й+і
Указанная методика позволяет решать задачу определения наиболее информативных показателей. Для этого столбцы матрицы А изначально следует заполнить данными различных показателей, располагая в каждом столбце данные одного срока наблюдения. Кроме того, в А можно поместить и всю информацию (по различным показателям и из разных пунктов наблюдения) так, чтобы номер столбца соответствовал общему временному отсчету, тогда можно выявить наиболее информативные данные с учетом вида показателя и места сбора информации.
Оценка информативности косвенных данных
Кроме непосредственных данных в экологическом мониторинге могут использоваться и косвенные источники информации, такие как, например, медицинские данные по заболеваемости населения и т.п. Для этих показателей (в общем количестве к)
можно построить другую матрицу В и, соответственно, В по той же методике (см. выше). Таким образом, количество строк в этой матрице равно к, а количество столбцов в будет равно п, как и в А.
Пусть вычислена матрица Z=P'h А размерности кхп и — BZ'А^ 1 размерности
п
тхк, тогда имеет место регрессионная модель х=¥1+е, где под вектором х подразумевается произвольный столбец В (набор центрированных и нормированных данных косвенной информации), а е — вектор погрешностей регрессии, обладающий теми же свойствами, что и е (см. выше). Далее, используя количество информации Is(xu,z)=У2logbdet(I+ЛhМ(u)), где ие{ 1,2,.,к}, М(и)=ае—^'(и^(и), мы можем так же, как и в предыдущем разделе, определить наиболее информативные косвенные показатели по отношению к общей экологической обстановке в анализируемом регионе (т.е. по отношению к факторам, ее определяющим). Смещенная оценка дисперсии компонент вектора погрешностей регрессии е в этом случае может быть получена по формуле
Ое2=(пк)~11т( В —FZ)( В —FZ)'
Описанный комплекс вычислительных процедур прошел успешную апробацию на данных медико-экологического мониторинга Архангельской области, собранных в 1996 г. за шестилетний период наиболее трудного этапа перестройки. Анализ проводился под эгидой Санкт-Петербургской государственной педиатрической медицинской академии, где на тот момент функционировала лаборатория, занимающаяся медико -экологическим мониторингом. Несмотря на хорошие с точки зрения вполне разумной интерпретируемости достаточно содержательные результаты анализа, руководство проекта отказалась от публикаций. Проведенная оценка информативности данных медицинского обеспечения (служивших косвенным источником информации) относительно данных по заболеваемости свидетельствовала о падении авторитета медицины у населения области, имея в виду амбулаторное обслуживание. Это отразилось в максимальной информативности данных скорой помощи относительно общей заболеваемости (те. заболевший человек сам не пытался обратиться к врачу) и лишало автора возможности указать опубликованный пример эффективного использования метода, рассмотренного в последнем разделе.
Общее замечание: Очевидно, что рассмотренная схема анализа при практическом использовании должна периодически повторяться (обновляться), так как ничего абсолютно постоянного в природе не существует, а человеческая активность не добавляет устойчивости. Однако это обновление (повтор) уже не представит каких-либо
существенных затруднений, так как предполагает использование того же программного комплекса, но с обновленными исходными данными.
Литература
1. Айвазян С. А., Бухштабер В. М., Енюков И. С., Мешалкин Л. Д. Прикладная статистика. Исследование зависимостей. - М.: Финансы и статистика, 1985. - 487 с.
2. Айвазян С. А., Бухштабер В. М., Енюков И. С., Мешалкин Л. Д. Прикладная статистика. Классификация и снижение размерности. - М.: Финансы и статистика, 1989.- 607 с.
3. Гельфанд И.М., Яглом А.М. О вычислении количества информации о случайной функции, содержащееся в другой такой функции // Успехи математических наук. 1957. Т. 12, № 1. - С. 3-52.
4. Пичугин Ю.А. География динамической неустойчивости циркуляции атмосферы в Северном полушарии (Моделирование и анализ) // Изв. РГО. 2005. Т. 137. Вып. 3. - C. 12-16.
5. Пичугин Ю.А. Замечания к отбору данных в задачах, связанных с регрессией // Заводская лаборатория. Диагностика материалов, 2002. № 5. - C. 61-62.
6. Себер Дж. Линейный регрессионный анализ. - М.: Мир, 1980. - 456 с.
7. Чанга Л.Б., Пичугин Ю.А. Распределение осадков и температуры воздуха в Танзании // Изв. РГО. 2006. Т. 138. Вып. 4. - C. 54-61.
8. Jackson D: Stopping rules in principal components analysis: a comparison of heuristical and statistical approaches. Ecology 1993, 74(8):2204-2214.
9. Bartkowiak A: How to reveal the dimensionality of the data? Applied Stochastic Models and Data Analysis 1991, 55-64.