Опыт применения геостатистического подхода к прогнозированию

Введение

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

В данной статье описана попытка применения геостатистических методов[2,3,4] для прогнозирования.

В качестве примера приводятся результаты ретроспективного прогноза максимальных уровней оз. Ильмень в период весеннего половодья, выполненные геостатистическим методом. Для сравнения приведены результаты детерминистского прогноза методом, описанным в [1] и, для нескольких лет, по методике ВФ ГГИ[7] 1994 года.

Провести эксперимент на данных по водному режиму оз. Ильмень предложил с.н.с. ВФ ГГИ А.В.Кокорев. Он также оказывал всемерную консультативную помощь и предоставил данные по режиму озера и его притоков.

Согласно А.П.Резникову[5], существуют "индуктивный и дедуктивный подходы к прогнозированию естественных процессов."

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

"При дедуктивном подходе специалисты в области методов обработки информации, руководствуясь довольно общими сведениями о свойствах целого класса процессов, создают методы прогнозирования, т.е. способы построения моделей прогнозирования, как правило, обширного класса процессов. Уже имея метод, при решении конкретной задачи наполняют модель конкретным содержанием, исходя из теоретических представлений о данной задаче и имеющихся фактических сведений. Таким образом двигаются от метода через модель конкретного процесса к самому акту прогнозирования"[5]. Прогноз строится не как функция факторов, а как функция наблюденных ранее значений предиктанта. Значения факторов используются при вычислении параметров этой функции.

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

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

Набор n значений факторов можно рассматривать как точку в n-мерном факторном пространстве. Каждой такой точке соответствует значение подлежащей прогнозу характеристики (в терминах метода - ПСФ). Приступая к задаче прогнозирования, надо располагать достаточно репрезентативным набором наблюдений, относящихся к различным моментам времени. Само время может рассматриваться как фактор, а может и не рассматриваться. В последнем случае "прогнозность" задачи оценивания ПСФ существует только в сознании прогнозиста, практически же выполняется интерполяция, или, реже, экстраполяция ПСФ в пространстве факторов, имеющих смысл физических причин. Если время входит в состав факторов, то при выполнении прогнозов всегда происходит экстраполяция. Далее будем считать, что время в состав факторов не входит и, если не будет необходимости особо подчеркнуть факт выхода точки прогнозирования за выпуклую оболочку исходных точек, ограничимся термином "интерполяция", понимая его в широкои смысле.

Теперь можно конкретизировать: под геостатистическим методом в данном тексте понимается геостатистический метод интерполяции.

Геостатистический метод можно рассматривать как альтернативу обучающейся прогнозирующей системе (ОС)[6].

Сильно упрощая, можно сказать, что ОС на этапе обучения решает задачу прямой оптимизации параметров прогнозной модели, для чего, многократно изменяя эти параметры и выполняя ретроспективные прогнозы, старается получще подогнать прогнозные оценки к факту.

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

Для первых опытов прогноза нами был выбран наиболее простой, т.н. ординарный (линейный точечный) кригинг, своиства которого хорошо известны.

В качестве анализа пространственной изменчивости ПСФ выполнялся вариограммный анализ, и далее будет упоминаться только он.

Геостатистический метод интерполяции

Геостатистика как совокупность математически обоснованных методов появилась в связи с задачами оценки содержания и запасов руд. Основы теории заложены Ж.Матероном[4]. В частности, им разработана теория пространственной случайной функции (ПСФ) и указаны ограничения на ее свойства, позволяющие обосновать применение разработанных расчетных методов. Отдельные методы применялись и до появления работ Матерона, в частности горным инженером Криге, в честь которого и назван основной расчетный метод геостатистики - кригинг.

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

Детальное изложение методов геостатистики выходит за рамки статьи. Их теоретические основы и расчетные схемы подробно описаны в литературе[2,3,4,9], что позволяет останавливаться лишь на моментах, важных для общего понимания.

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

Линейный кригинг вычисляет оценку ПСФ в указанной точке как линейную комбинацию значений ПСФ в исходных точках. Предварительно решается т.н. система уравнений кригинга относительно коэффициентов линейной комбинации.

Множество кригинг-оценок ПСФ в точках пространства образуют интерполирующую гиперповерхность(ИГ). Будем различать строгую интерполяцию и сглаживающую (аппроксимацию). При строгой интерполяции значения ПСФ в исходных точках принадлежат ИГ. При аппроксимации ИГ лишь примерно описывает ПСФ в исходных точках, за счет чего можно достичь дополнительной гладкости ИГ.

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

Точечный кригинг - строго интерполяционная процедура. Матрица системы уравнений кригинга при совпадении исходных точек вырождается.

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

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

Вариограмма определяет форму ИГ в промежутках между исходными точками. Влияние вариограммы на ИГ больше в тех областях пространства, где исходные точки разрежены и наиболее сильно в области экстраполяции.

Вариограмма также дает представление о таком важном свойстве ПСФ, как анизотропия.

Численное значение вариограммы при нулевом векторе приращений координат в геостатистике называют эффектом самородков(ЭС). Его значение лежит в закрытом интервале [0, D], где D - общая дисперсия значений ПСФ. Большее нуля значение ЭС говорит о неоднозначности ПСФ.

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

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

Специфика и проблемы

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

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

С ростом числа факторов увеличивается и потребный объем исходных данных, что вынуждает применять локальную интерполяцию во избежание чрезмерных порядков системы уравнений кригинга. В связи с этим возникает проблема метода поиска опорных точек, поскольку из существующих методов, избегающих экстраполяции, не пригоден, по большому счету, ни один. Это связано с тем, что, избегая экстраполяции, методы поиска набирают опорные точки с избытком, растущим с ростом размерности. При размерности факторного пространства не более 5-6-ти еще может применяться квадрантный метод. При использовании этого, наиболее экономного из известных, метода число опорных точек одного ореола составляет 2ⁿ. Для сравнения, минимальный ореол, страхующий от экстраполяции (симплекс) содержит n+1 точкy (n - размерность пространства). В связи с изложенным был разработан алгоритм последователного квазиоптимального включения опорных точек с учетом ковариации включаемой точки как с точкой оценивания, так и с включенными ранее.

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

А.П.Резников[6], считает аппроксимативность необходимым свойством прогнозирующей системы с точки зрения устойчивости к ошибкам в исходных данных.

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

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

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

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

Озеро Ильмень как объект прогнозирования

Озеро Ильмень расположено в Новгородской области, в наиболее пониженной части Приильменской низины на высоте около 18м над уровнем моря и является одним из крупнейших в Европе.

Площадь бассейна оз.Ильмень - 67300 кв. км. Большая часть ее относится к водосборам четырех главных притоков: рек Мсты, Ловати, Шелони и Полы. Сток из озера происходит в р.Волхов, на которой в 1928г. построена плотина Волховской ГЭС.

Минимальный уровень озера - 16.5м (1972 г.), при этом площадь озера составляет 750 кв. км. В половодье уровень может достигать 22.7м (1966 г.), что "приводит к затоплению пологих склонов Ильменской котловины на огромных пространствах"[7] - площадь озера увеличивается до 2100 кв.км, т.е в 2.8 раза. Объем призмы регулирования достигает 9 куб. км. (Данные приведены за период 1956 - 1990 гг.)

Ясно, что прогноз уровней оз.Ильмень в периоды половодья является не только интересной, но и важной задачей. В разное время для ее решения разрабатывались различные методики. Отметим методику бюро гидрологических прогнозов СЗ УГУКС, разработанную под руководством Р.Е.Нежиховского и наиболее современную методику ВФ ГГИ [7] 1994 года. Последняя строилась с учетом закрытия многих гидрометрических постов в связи с недостатком финансирования.

Остановимся более подробно на методике ВФ ГГИ, поскольку рассматриваемые нами геостатистические прогнозы выдаются в те же сроки с использование части той же информации.

Методика ВФ ГГИ предполагает выдачу прогноза через 1-2 дня после прохождения пика расходов на главных притоках озера. Прогнозируется ход уровня с пентадным шагом по времени на весь период половодья.

Максимум уровня на озере наступает, как правило, через 1-3 недели после прохождения пика паводка на притоках. Повышение уровня с момента выдачи прогноза до максимального в разные годы меняется почти от нуля до двух с лишним метров, чему соответствует увеличение площади озера на 500 кв. км.

Методика ВФ ГГИ предполагает поступление информации об уровнях оз.Ильмень и расходах 4-х главных притоков, посты на которых контролируют 75% площади бассейна озера. Для моделирования подпора со стороны р.Волхов используется информация о расходах притоков р.Волхов - рек Тигоды, Керести и Пчевжи, а также информация об уровне верхнего бьефа Волховской ГЭС. Могут привлекаться прогнозы осадков.

Для оценки параметров прогнозной модели ВФ ГГИ и ее верификации использовались данные за период 1956-1990гг.

Для проверки применимости геостатистического подхода и его программной реализации в "Геолинкеquot; использовался тот же фактический материал, но число информационных гидрометрических постов было уменьшено за счет постов на притоках р.Волхов и Волховской ГЭС и не использовались прогнозы осадков.

Были выполнены ретроспективные прогнозы, идентичные по форме прогнозам ВФ ГГИ и, отдельно, прогнозы максимальных уровней с тем же сроком выдачи.

Верификация прогнозов

Все рассматриваемые в данном разделе прогнозы, также как в методике ВФ ГГИ, выдаются после прохождения пика расхода основных притоков озера в створах р.Мста-Девкино, р.Ловать-Холм, р.Шелонь-Заполье и р.Пола-Налючи.

Было выполнено сравнение результатов геостатистического и детерминистского прогноза, построенного на основе уравнения водного баланса в предположении линейной зависимости прирашений уровня и объема, а также стока из озера и уровня[1].

В обоих случаях в качестве факторов использовались максимальный суммарный расход притоков Q_max и уровень Н₀ озера на момент максимума расхода.

Прогнозировались максимальные приращения уровней H_max-H₀, по которым легко вычислить прогноз значения H_max.

Использовались данные за период 1956-1990гг. В 1976, 1981 и 1990 годах рассматривались два максимума.

Среднеквадратические ошибки (СКО) обоих прогнозов примерно равны (16см и 17см). Корреляция прогнозных и фактических значений весьма высокая (r>0.9). Эти оценки соответствуют длине ряда наблюдений за паводками, равной 37. Прогнозы обоих типов ближе друг к другу, чем каждый из них к факту (см. рис.1 и рис.2). Можно полагать, что основным источником ошибок прогноза в обоих случаях является неучет фактора осадков в период между моментом выдачи прогноза и наступлением максимума уровня.

Для оценки точности прогнозирования обоими методами применялась процедура, известная в геостатистике под названием перекрестной проверки (cross validation [9]).

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

На рис.1 изображены графики прогнозов (точнее - перекрестных оценок) обоими методами и график фактических значений максимальных уровней.

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

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

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

На рис.3 представлен график зависимости СКО геостатистического прогноза от длины ряда наблюдений.

Рис.3 Зависимость СКО геостатистического прогноза от длины ряда наблюдений. Аппроксимация проведена от руки. Можно полагать, что СКО спадает по экспоненте с асимптотой на уровне квадратного корня из ЭС.

Чтобы выявить разницу в "поведении" прогнозных моделей, ослабили один из факторов: вместо суммарного расхода использовали только расход р.Мста в створе Девкино, контролирующего 30% площади бассейна. На рис.4 хорошо видно, что геостатистическая модель, опирающаяся только на данные (без поддержки разумной теорией), при малой длине ряда совершенно беспомощна. Она выходит на режим позже, чем модель на основе уравнения водного баланса и по мере накопления данных догоняет последнюю по качеству прогноза.

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

В отчете ВФ ГГИ приведены результаты прогноза уровня оз.Ильмень в периоды весеннего половодья 1984-1988 годов. Для этих лет в таблице 1 приведены ошибки прогноза максимальных уровней тремя методами. В данном случае прогнозы так же строились "честно", с использованием исключительно предшествующих наблюдений. Заметим, что по методике ВФ ГГИ прогнозировались не максимумы, а временные ряды, и значения H_max выбиралось как максимальное значение в ряду.

Сотрудницей "Геолинка" М.А.Эрастовой в ходе работы над дипломом с применением 1-ой версии программы были выполнены геостатистические ретроспективные прогнозы, идентичные по форме прогнозам ВФ ГГИ. Прогнозировался уровень озера с упреждением от 5 суток с пентадным шагом. Использовались три фактора - два указанных выше и дополнительные, различные для разных упреждений. Оценка точности этих прогнозов приведена в таблице 2, заимствованной из диплома[8]. В отчете ВФ ГГИ [7] статистическая оценка точности не приводится.

Заключение

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

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

К сожалению, анизотропный вариограммный анализ в многомерном пространстве вызывает наибольшие затруднения. Трудозатраты на программную реализацию, видимо, сравнимы с таковыми при разработке обучающейся системы(ОС). Надежды на существенную экономию трудозатрат при выборе геостатистического подхода как альтернативы ОС не оправдались.

Отметим следующий обидный парадокс. С ростом длины ряда наблюдений моделирование вариограммы облегчается и растет надежность оценок ее параметров. Одновременно растет густота заполнения области изменения факторов точками наблюдений и уменьшается влияние деталей формы аппроксимирующей поверхности, определяемых вариограммой, на результат. Изощренный аппарат вариограммного анализа с его многомерными оптимизациями становится все менее нужным. Мы убедились в этом при построении прогнозов суточной аккумуляции воды озером с использованием ряда наблюдений длиной около 13000.

Разработанный в "Геолинке" программный аппарат в настоящее время еще нельзя считать законченым продуктом, работа продолжается.