«В. В. Демьянов, Е. А. Савельева ГЕОСТАТИСТИКА теория и практика Под редакцией профессора, доктора физико-математических наук Р. В. Арутюняна Москва Наука 2010 УДК 91:519.8 ББК 26.8в6 Г35 Рецензенты: доктор технических ...»

7. Многоточечное моделирование (multipoint statistics simulation) — первоначально реализовано для категориальных данных, но недавно были разработаны алгоритмы многоточечной статистики для моделирования и непрерывных переменных. Этот подход подробно рассмотрен в Главе 11.

Стохастическое моделирование пространственной неопределенности 8.2. Последовательный принцип моделирования Геостатистика интерпретирует измерения пространственной переменной Z(xi) в точках xi (i = 1,..., n) как реализацию значений случайной функции Z(x), которая определенна в области S и характеризуется совместной условной функцией распределения F(x|z). В рамках такой системы проблема заключается в генерации K реализаций этой случайной функции, детально покрывающих область S набором из N точек [x1,..., xN]. Каждая реализация соответствует совместной условной функции распределения, в которую включены все статистические характеристики процесса:

Принцип последовательного моделирования основан на использовании правила Байеса:

Применяя последовательно (N раз) формулу (8.2) к правой части формулы (8.1) и переписав ее в виде функции распределения, получаем:

где F ( x N 1 ; z N 1 | ( n + N 2) ) — условная функция распределения Z(x), определяемая n исходными данными и N – 1 значениями реализации z(xj), j = 1,..., N – 1. Полученное выражение и является теоретической основой последовательных методов стохастического моделирования.

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

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

1. Построение случайной последовательности из набора N [x1,..., xN] точек для оценки x1,..., x N, где k — номер реализации.

2. Для точки x1 проводится оценка локальной функции распределения по набору исходных данных. В соответствии с оцененной функцией расВ. В. Демьянов, Е. А. Савельева Геостатистика: теория и практика пределения разыгрывается значение функции Z(x) в точке x1k z x1k.

Это значение добавляется к исходным данным.

3. Для точки xm проводится оценка локальной функции распределения по набору исходных данных и полученным на предыдущих m – 1 шагах значениям z x1k,..., z xm 1. В соответствии с оцененной функцией расk пределения разыгрывается значение функции Z(x) в точке xm z xm. k Это значение также добавляется к исходным данным.

4. Шаг 3 повторяется для всех последующих N – 1 точек в соответствии с последовательностью, полученной на шаге 1. Каждый раз оценка локальной функции распределения производится по исходным данным и значениям, сгенерированным на предыдущих шагах.

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

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

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

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

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

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

Многие модели геостатистического последовательного стохастического моделирования используют кригинг при моделировании локальной функции распределения: либо для моделирования параметров распределения, либо для моделирования табулированных вероятностей для непараметрических методов. При использовании кригинга может возникнуть эффект экранирования [Isaaks, Srivastava, 1989]. Он состоит в уменьшении веса точек, попадающих между одной из точек измерения и точкой оценивания, что может привести к появлению отрицательных весов. Таким образом, эффект экранирования будет усиливаться за счет добавления вновь смоделированных значений к данным после каждого шага оценивания. На практике не обязательно использовать все существующие значения измерений для построения условного распределения в точке оценивания. Необходимо ограничить используемые условные измерения окрестностью точки оценивания, при этом можно ограничивать отдельно количество исходных данных и количество уже смоделированных значений. Более подробно эффект экранирования описан в Главе 5.

Среди наиболее широко используемых алгоритмов последовательного моделирования можно выделить:

• гауссово;

• обрезанное гауссово;

• индикаторное;

• прямое;

• многоточечное и др.

Первые два алгоритма являются параметрическими, остальные — непараметрическими. Эти алгоритмы отличаются различными предположениями о локальном законе распределения.

Далее мы рассмотрим некоторые наиболее часто используемые алгоритмы.

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

В этом случае для любой точки области локальная функция распределения будет распределена по нормальному закону и будет определяться двумя параметрами — средним и вариацией.

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

На первом этапе моделирования предполагается стационарность случайной функции Z(x) и существование такой случайной функции Y(x) где — однозначная функция нормализующего (normal score) преобразования (рис. 8.3).

В. В. Демьянов, Е. А. Савельева Геостатистика: теория и практика Гауссово преобразование производится путем постановки в соответствие реальной кумулятивной функции распределения исходных данных кумулятивного стандартного нормального распределения:

где F() — кумулятивное распределение частоты исходных данных; G() — стандартное кумулятивное нормальное распределение. Такое соответствие проиллюстрировано на рис. 8.3. На практике при проведении прямого преобразования строится специальная таблица соответствия значений для обратного преобразования.

После преобразования распределение данных становится нормальным, т. е.

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

где C(h) — ковариация; (h) — вариограмма функции Y(x).

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

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

Стохастическое моделирование пространственной неопределенности Для проверки на бинормальность условной функции распределения любого набора пар данных {y(xi), y(xi + h), i = 1,..., N(h)} используется ковариация CY(h). Существуют аналитическое и табулированное соотношения между ковариацией CY(h) и значением стандартной нормальной функции распределения [Deutsch, Journel, 1998]:

где yp = G–1(p) — стандартный нормальный p-квантиль; CY(h) — ковариация стандартной нормальной случайной функции Y(x).

Бинормальная вероятность является нецентральной индикаторной ковариацией (см. (7.3)) для порога yp:

где I(x, p) — индикаторное преобразование для функции Y(x) (см. (7.1));

I(h, p) — индикаторная полувариограмма для p-квантиля отсечения yp.

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

Существует и более грубый тест на бинормальность [Emerly, 2005]. Он заключается в проверке соотношения вариограммы и мадограммы М:

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

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

В. В. Демьянов, Е. А. Савельева Геостатистика: теория и практика Рис. 8.4. Проверка на бинормальность исходных данных Следующим шагом является построение ковариации для нормализованных значений Y(x) и ее моделирование. Эта задача подробно описана в Главе 4.

Напомним, что ковариация переменной y для вектора h = xi – xj обозначается как Cij = C(xi – xj).

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

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

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

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

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

3. В каждой точке оценивания производятся:

Стохастическое моделирование пространственной неопределенности • оценка параметров локального распределения плотности вероятности (среднего и вариации) с помощью простого кригинга на основе преобразованных исходных данных и уже сгенерированных значений в других точках последовательности;

• выборка случайного значения в соответствии с нормальным распределением и полученными параметрами;

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

4. Обратное преобразование промоделированных реализаций из нормализованных значений.

Для получения следующей равновероятной реализации повторяются шаги 2—4.

Нормализация исходных данных Использование смоделированного Обратное преобразование нормализованных значений Рис. 8.5. Схема алгоритма последовательного гауссова моделирования В. В. Демьянов, Е. А. Савельева Геостатистика: теория и практика Ключевой момент последовательного гауссового моделирования — построение в каждой точке оценивания нормальной функции плотности вероятности. Локальная гауссова функция распределения определяется двумя параметрами — средним значением и вариацией. Для их получения используется простой кригинг с известным средним значением. Обычный кригинг в движущемся окне может быть использован для оценки среднего значения в нестационарном случае, но при этом вариация оценивается простым кригингом.

Оценка простого кригинга YSK ( x ) в точке x вычисляется так (см. Раздел 5.2):

где среднее значение нормализованных данных m по области постоянно (в предположении о стационарности второго порядка) и равно нулю.

Веса простого кригинга определяются путем решения системы n(x) уравнений простого кригинга:

где n(x) — общее количество использующихся соседних точек xi при оценке точки x.

Вариация оценки простого кригинга В результате получаем параметры локальной нормальной функции распределения N YSK ( x ), 2 ( x ) в точке оценивания x.

Далее проводится случайная выборка (по методу Монте-Карло) из полученного нормального распределения. Разыгранное значение является равновероятной стохастической реализацией значения функции Y(x) в данной точке.

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

После прохода через все точки оценивания для получения окончательного результата моделирования полученные нормальные значения оценок { y ( x ), x A} преобразуются обратно в абсолютные значения исходной Стохастическое моделирование пространственной неопределенности функции z ( x ) = 1 ( y ( x ) ), x A с использованием обратного гауссова преобразования.

Если модель хорошо соответствует исходным данным, то реализации последовательного гауссова моделирования воспроизводят:

• стандартное нормальное распределение у преобразованной переменой;

• вариограммы преобразованных переменных;

• значения преобразованных переменных в точках измерений.

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

Приведенная выше выборка из нормального распределения со средним, равным оценке простого кригинга YSK ( x ), эквивалентна случайной выборке из нормального распределения с нулевым средним N ( 0, SK ( x ) ).

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

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

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

В. В. Демьянов, Е. А. Савельева Геостатистика: теория и практика Примеры реализаций последовательного гауссова моделирования с различными параметрами вариограммы приведены на рис. 8.6 для уровня наггета 0 и 40% от априорной вариации (а, б); изотропии (г, е), геометрической анизотропии (в) и зонной анизотропии (д) радиуса корреляции; угла направления длинной корреляционной структуры от вертикали 0° (ж), 30° (в) и 90° (з); соотношения радиусов корреляции в ортогональных направлениях 8 и 40 (в), 80 и 4 (д), 40 и 8 (ж, з) в единицах расстояния.

Стохастическое моделирование пространственной неопределенности Рис. 8.6. Примеры стохастических реализаций последовательного гауссова моделирования для различных значений параметров модели вариограммы:

а — R = 40, r = 8, наггет = 0.0, угол = 60; б — R = 40, r = 8, наггет = 0.4, угол = 60;

в — R = 8, r = 40, наггет = 0.0, угол = 60; г — R = 40, r = 40, наггет = 0.0, угол = 60;

д — R = 80, r = 4, наггет = 0.0, угол = 60; е — R = 8, r = 8, наггет = 0.0, угол = 60;

ж — R = 40, r = 8, наггет = 0.0, угол = 0; з — R = 40, r = 8, наггет = 0.0, угол = Предположение о локальной нормальности имеет и негативные стороны. Гауссово моделирование является алгоритмом максимальной энтропии — максимального беспорядка в стохастической реализации. Это означает слабую связанность предельных значений, т. е. точки с максимальными значениями переменной не будут иметь связи друг с другом по соседним ячейкам с высокими значениями переменной. Такое поведение не характерно, например для В. В. Демьянов, Е. А. Савельева Геостатистика: теория и практика геологических приложений, где пласты высокой проницаемости образуют связные структуры. Дело в том, что вариограмма характеризует корреляцию, основываясь на связи пар точек, в то время как связная структура, образованная несколькими точками, не сохраняется и может быть разрушена при стохастическом моделировании, что и происходит в случае максимальной энтропии (беспорядка), свойственной гауссову моделированию.

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

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

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

Локальное нормальное распределение непрерывной переменной, полученной в результате последовательного гауссова моделирования, можно разбить по категориям на основе выбранных пороговых значений (рис. 8.7). Так определяются значения искомой категориальной (или разрывной) переменной. Значение стохастической реализации получается при попадании случайно выбранного значения из локального нормального распределения в тот или иной интервал.

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

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

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

Индикаторный подход заключается в моделировании бинарных индикаторных переменных, которые принимают значения либо 1, либо 0 в зависимости от присутствия или отсутствия свойства в данной точке. Индикаторный подход может быть использован для моделирования как категориальных, так и непрерывных переменных. Для получения значений бинарных индикаторных переменных проводится индикаторное преобразование исходных данных для выбранного набора срезов в случае непрерывной функции (7.1) В случае категориальной переменной значения индикаторных переменных (1 и 0) соответствуют присутствию или отсутствию каждой из категорий в В. В. Демьянов, Е. А. Савельева Геостатистика: теория и практика точке измерения — см. (7.2). Индикаторные преобразования, выбор числа и значения срезов обсуждались в Разделе 7.1.

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

где pk = E I ( x; sk ) [0, 1] определяет долю категории zk в глобальном распределении и находится из исходных данных с учетом декластеризации.

Веса i определяются индикаторным кригингом с использованием модели ковариации для соответствующих категорий zk (см. Главу 7). Если средние значения доли категории варьируются по области, то можно использовать простой индикаторный кригинг с гладко меняющимся локальным средним значением.

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

Алгоритм последовательного индикаторного моделирования заключается в следующих этапах [Goovaerts, 1997] (рис. 8.8).

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

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

3. В каждой точке оценивания моделируется стохастическая реализация по следующей последовательности операций (см. рис. 8.7):

• оценка K вероятностей Pk(x|(zk)) k = 1, …, K при помощи индикаторного кригинга в выбранной точке x последовательности;

• построение локальной условной функции плотности вероятности на основе K вероятностей Pk (коррекция, интерполяция, экстраполяция, как описано в Главе 7);

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

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

4. Переход к следующей точке оценивания и выбранной последовательности и повторение шагов 2 и 3.

Шаги 2—4 повторяются для получения нескольких равновероятных реализаций в точках оценивания.

преобразование Случайный выбор ячейки сетки оценивания.

Значения двух бинарных индикаторных переменных Оценка индикаторного кригинга вероятности: категориального значения либо превышения для каждого порогового значения Построение локальной функции плотности распределения по вероятностям, оцененным индикаторным кригингом, и случайная выборка промоделированного значения Использование смоделированных значений индикаторных переменных при оценке последующих ячеек Рис. 8.8. Схема алгоритма последовательного индикаторного моделирования В. В. Демьянов, Е. А. Савельева Геостатистика: теория и практика Ни рис. 8.8 изображена пошаговая схема индикаторного моделирования на примере двух категориальных переменных, которым соответствуют индикаторные переменные. Значения переменных в точках даны вне скобок и в скобках.

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

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

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

Стохастическое моделирование пространственной неопределенности Рис. 8.9. Реализации последовательного индикаторного моделирования для различных горизонтальных радиусов корреляции:

(вертикальный радиус корреляции R = 8) В. В. Демьянов, Е. А. Савельева Геостатистика: теория и практика Рис. 8.10. Реализации последовательного индикаторного моделирования для различных вертикальных радиусов корреляции:

Индикаторное моделирование отличается следующими особенностями.

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

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

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

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

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

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

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

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

В. В. Демьянов, Е. А. Савельева Геостатистика: теория и практика 8.6. Последовательное прямое моделирование При прямом моделировании отсутствует предварительное преобразование исходной переменной с последующим поиском функции распределения в преобразованном пространстве, что характерно для двух описанных выше методов.

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

Итак, мы рассматриваем непрерывную функцию Z(x) с глобальной условной функцией распределения FZ ( z ) = Pr Z ( x ) < z и стационарной вариограммой (h), которые мы заинтересованы воспроизводить в наших реализациях. При этом для воспроизведения вариограммы нам достаточно использовать локальные условные кумулятивные функции распределения плотности вероятности, центрированные в оценке простого кригинга:

где xi — местоположения данных (исходных и уже смоделированных) с условными вариациями, полученными из простого кригинга 2 ( xu ). При этом неважно, какая именно функция распределения используется.

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

Стохастическое моделирование пространственной неопределенности Более удобный подход состоит в использовании разбиения на классы, аналогичного обрезанному гауссову распределению [Deutsch, 2002].

Преобразуем исходные данные z(x) в нормальное распределение при помощи функции :

Локальная оценка простого кригинга z*(xu) имеет эквивалентное гауссово значение y*(xu) = (z*(xu)), которое совместно со стандартизованной вариацией простого кригинга SK ( xu ) позволяет определить гауссову функцию распределения G y ( xu ), 2 ( xu ).

Эта гауссова условная функция распределения позволяет определить интервал условной функции распределения z(x), в котором нужно разыгрывать новое значение:

• сгенерировать значение p из распределения U(1, 0);

• сгенерировать значение ys из распределения G y * ( xu ), 2 ( xu ) :

Разыгрываемое значение получается обратным преобразованием:

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

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

Примеры, представленные на рисунках 8.11, 8.12, иллюстрируют это качество метода. Рассмотрены исходные данные с двумя типами глобальной функции распределения: характеризующейся двумя пиками (см. рис. 8.11) и примерно однородным (см. рис. 8.12). Реализации, различающиеся между собой, достаточно четко воспроизводят исходные гистограммы (рис. 8.13 — 2 пика, рис. 8.14 — однородное распределение).

В. В. Демьянов, Е. А. Савельева Геостатистика: теория и практика Рис. 8.11. Исходные данные и их глобальная функция распределения Рис. 8.12. Примеры реализации прямого условного моделирования Стохастическое моделирование пространственной неопределенности Рис. 8.13. Исходные данные и их глобальная функция распределения Рис. 8.14. Примеры реализации прямого условного моделирования и их глобальные функции распределения В. В. Демьянов, Е. А. Савельева Геостатистика: теория и практика 8.7. Моделирование отжига Моделирование отжига (simulated annealing) [Metropolis et al., 1985] — это общее название для семейства оптимизационных алгоритмов, основанных на принципах стохастической релаксации. Моделирования отжига аналогично процессу остывания металла в термодинамике и основано на соотношении Больцмана между температурой и энергией [Aarts, Korst, 1989].

Термин «отжиг» (annealing) пришел в математику из металлургии. При высокой температуре металл легко деформируется и меняет форму. Чем выше температура, тем больше скорость колебаний атомов и тем легче металл поддается деформации. Если резко снизить температуру, то атомы «замрут»

и мы получим твердый, но очень хрупкий металл. Если стоит задача получить максимально организованную структуру, необходимо сначала сильно разогреть металл, а затем очень медленно охлаждать его. Таким образом, металл будет проходить через множество квазиравновесных состояний и у каждого его атома будет достаточно времени, чтобы найти «лучшее» место среди других атомов в смысле минимума полной энергии системы, что соответствует идеальной кристаллической решетке. С точки зрения математики задачу можно сформулировать так: минимизировать среднюю квадратичную ошибку отклонения уровня поверхности металла от некоторого постоянного значения. Оптимизируемыми параметрами в этом случае будут координаты положения атомов.

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

Рассмотрим моделирование непрерывной величины z в N узлах сетки xj при заданных условиях z(x), = 1,..., n таким образом, чтобы вариограмма данных воспроизводилась для первых S лагов. Аннилинг требует задания целевой функции (являющейся аналогом энергии), которая измеряет разницу между значениями целевых и текущих статистических параметров на каждом i-м возмущении. Если цель — воспроизвести вариограммную модель, то целевая функция может выглядеть так:

Стохастическое моделирование пространственной неопределенности где (hs) — значение требуемой вариограммной модели для лага hs и (i (hs) — соответствующее значение вариограммы реализации на i-м возмущении.

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

Общий алгоритм моделирования отжига состоит из следующих этапов.

1. Создаем начальный образ z( 0 ) x j, j = 1,..., N, который сохраняет исходные данные и может сразу аппроксимировать какую-нибудь целевую статистику (дисперсию распределения, плато вариограммы или гистограмму).

2. Считаем начальное значение целевой функции, соответствующее этой начальной реализации.

3. Возмущаем реализацию каким-либо механизмом, например отражением пар z-значений: z(0)(xj) становится z(0)(xi), и наоборот. По аналогии с физическими процессами в остывающем металле механизм возмущения зависит от температуры: чем она ниже, т. е. чем больше шагов сделано, тем меньше меняется значение в точке при возмущении. Например, в случае отражения пар уменьшается расстояние между точками, значения которых мы меняем местами.

4. Оцениваем эффект возмущения на воспроизведение целевой функции, снова вычисляя ее значение Onew, учитывая модификацию начальной реализации.

5. Принимаем или не принимаем возмущение на основе какого-либо правила. Обычно вероятность принятия задается с помощью распределения Больцмана [Aarts, Korst, 1989]:

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

6. Если возмущение принимается, заменяем начальную реализацию на новую z(l1) x j, j = 1,..., N с соответствующей целевой функцией Оold = Onew.

В. В. Демьянов, Е. А. Савельева Геостатистика: теория и практика 7. Повторяем шаги с 3-го по 6-й, пока целевая структура не будет приемлемо достигнута или пока возмущения не перестанут уменьшать целевую функцию. Потом снижаем температуру и проделываем указанную процедуру (шаги 3—6) для новой температуры и так поступаем до тех пор, пока не достигнем приемлемого результата (см. ниже критерий остановки).

Последующие равновероятные реализации z ( l ) x j, j = 1,..., N, l l производятся повторением шагов 1—8, начиная с другого начального образа.

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

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

Требования к начальному образу.

• Он должен быть легко производим.

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

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

Обычно начальный образ производится так: исходные данные «замораживаются» на своих местах, а приписываемое каждому узлу значение z-величины выбирается случайным образом в соответствии с глобальной функцией распределения F(z). Такой подход достаточно быстр и дает набор начальных образов, уже удовлетворяющих целевой гистограмме.

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

Стохастическое моделирование пространственной неопределенности Чаще всего используются два механизма возмущения.

Первый — отражение z-величин для случайным образом выбранных пар точек xj и xk, находящихся на расстоянии d D(T), D(T), уменьшается понижением Т:

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

Другой механизм называют возмущением. Случайным образом выбирается одна точка xj, и модифицируется соответствующее значение z(i–1)(xj) согласно какому-нибудь механизму, например где pj — случайное значение из {0, 1}; а F(z) — целевая гистограмма.

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

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

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

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

Кумулятивная функция распределения (учет гистограммы данных).

Типичная целевая статистика — это по возможности декластеризованная однопеременная функция распределения z-данных F(z). Если диапазон изВ. В. Демьянов, Е. А. Савельева Геостатистика: теория и практика менения z описывается серией из К порогов zk, то разница между целевой и текущей функцией распределения может быть измерена так:

где F() ( zk ) — значение на пороге zk, вычисленное для реализации на i-м возмущении.

Модель полувариограммы. Воспроизведение вариограмной модели (h) обычно ограничивается определенным количеством S лагов. Разница между целью и текущим значением в этом случае измеряется так:

где g ( i ) (hs ) — значение полувариограммы реализации для лага hs на i-м возмущении.

Модели индикаторных полувариограмм. Моделирование отжига позволяет учитывать специальные пространственные структуры, моделируемые индикаторными полувариограммами I(h, zk), посчитанными для К различных порогов zk. Так получаем еще одну компоненту целевой функции:

где g (Ii ) ( hs ; zk ) — значения индикаторной полувариограммы на лаге hs по порогу zk для реализации на i-м возмущении.

Коэффициент корреляции. Пусть Y — лучшая выборка или ранее промоделированная величина, вторичная по отношению к величине Z. Если взаимосвязь Z-Y корректно описывается линейным коэффициентом корреляции rZY(0), то его можно ввести в целевую функцию в виде компоненты где (ZY (0) — коэффициент корреляции, вычисленный по соответствующим парам y-данных и z-значений на i-м возмущении.

Кросс-вариограмма. Пространственная кросс-корреляция между величинами Z и Y, моделируемая с помощью кросс-вариограммы ZY(h), может быть воспроизведена включением в целевую функцию компоненты типа Стохастическое моделирование пространственной неопределенности где g (ZY) ( hs ) — значение кросс-вариограммы между z-значениями и y-данными на i-м лаге.

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

• целевая функция достигла достаточно малого значения Omin;

• количество возмущений при одной температуре превысило допустимое максимальное число;

• доля приемлемых возмущений меньше, чем заданное пороговое значение.

Примеры реализаций, полученных с использованием моделирования отжига, приведены на рис. 8.15. Это реализации пространственного распределения краба Берди в 2006 г. в Беринговом море. Использование моделирования отжига в данном случае обусловлено необходимостью использования индикаторных вариограмм как более робастных (это показало в Главе 7).

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

Рис. 8.15. Примеры реализаций с использованием моделирования отжига для данных траловой съемки по пространственному распределению краба Берди В. В. Демьянов, Е. А. Савельева Геостатистика: теория и практика 8.8. Объектное моделирование Объектное моделирование является альтернативой пиксельному моделированию. В отличие от многих геостатистических моделей оно не основано на двухточечной статистике (вариографии). Однако объектный подход позволяет промоделировать пространственную корреляцию без помощи вариограммы на основе набора пространственных структур, которые имеют определенную заданную форму, — объектов. Для распределения объектов используются булево моделирование и алгоритмы оптимизации. Как правило, объектное моделирование применяется для категориальных переменных. При этом распределение рассматриваемой категории моделируется как совокупность геометрических объектов, которые покрывают области, в которых преобладают значения данной категории.

Примеры объектного моделирования Для проведения объектного моделирования требуется определить набор форм объектов. Обычно при этом руководствуются экспертным анализом на основе физических представлений об исследуемой системе. Так, в геологических приложениях формы объектов могут быть получены в результате анализа выходов пород, разрезов, шурфов, информации из буровых скважин, а также геологических представлений. Примеры геологических объектов — флювиальные синусоидные и меандрированные русла, песочные линзы, эоловые дюны нетривиальной трехмерной геометрии, диски или эллипсоиды, сланцевые останцы, конусы прорыва, лагуны и т. п. (рис. 8.16).

Рис. 8.16. Примеры объектного моделирования:

а — система речных русел с двумя типами объектов: связанные кривые русла с прирусловым валом и дугообразные формы старых русел;

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

Объекты могут размещаться в пространстве при помощи различных алгоритмов [Deutsch, 2002].

Наиболее простым методом является случайное размещение объектов.

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

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

При этом следует иметь в виду, что если точки измерения выбирались с предпочтением к определенным значениям переменных, то доля рассматриваемой категории в распределении исходных данных измерений может быть завышена (рис. 8.17) Расположение условных данных измерений Расположение объектов случайным образом в окрестности точек измерения так, чтобы они удовлетворяли данным измерений Расположение дополнительных объектов в областях отсутствия данных для достижения заданной доли значений категориальной переменной Рис. 8.17. Алгоритм объектного булева моделирования В. В. Демьянов, Е. А. Савельева Геостатистика: теория и практика Оптимальным и наиболее трудоемким методом является итерирование перемещений и изменение размеров объектов для наилучшего удовлетворения данных измерений и заданной доли категории. Для этого осуществляют минимизацию целевой функции при помощи стандартных алгоритмов (градиентных и стохастических), используя случайные возмущения в размещении объектов.

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

При использовании слишком большого разнообразия форм и степеней свободы объектов объектная модель может оказаться слишком громоздкой.

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

Еще одним недостатком объектного подхода является сложность включения в него дополнительных вероятностных «мягких» данных [Caers, 2005].

«Мягкие» (soft) данные обычно имеют вид вероятности значения переменной, которая определяется на регулярной сетке на основе данных более низкого разрешения, чем сетка оценивания. Эти данные коррелированны с моделируемой переменной, но требуют калибровки и менее точны (например, сейсмическое зондирование, аэрогаммасъемка). В случае пиксельных моделей такие данные приводятся на сетку оценивания более высокого разрешения с учетом изменения вероятности. При моделировании объектами, имеющими достаточно низкое разрешение, использование таких «мягких»

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

Упражнение 8.1. Почему оценку кригинга называют сглаженной, чем от нее отличается реализация стохастического моделирования?

Упражнение 8.2. При оценивании часто точно неизвестен верхний предел оцениваемой переменной. Чему равна максимальная/минимальная оценка кригинга? Что больше — максимальная оценка кригинга или максимальное значение стохастической реализации?

Стохастическое моделирование пространственной неопределенности Упражнение 8.3. Вариограмма характеризует пространственную корреляцию и уровень пространственной вариабельности. На рис. 8.18 приведены две вариограммы, построенные для оценки кригинга и стохастической реализации на основе одних и тех же данных и модели вариогрммы. Какая из вариограмм соответствует оценке кригинга, а какая — реализации стохастического моделирования? На чем основан выбор?

Рис. 8.18. Вариограммы для оценки кригинга и стохастической реализации Упражнение 8.4. В последовательном гауссовом моделировании используется модель вариограммы. Какие типы моделей вариограмм могут быть использованы? Почему суммарное плато модели вариограмм для гауссова моделирования должно быть равно единице?

В. В. Демьянов, Е. А. Савельева Геостатистика: теория и практика Литература Каневский М. Ф., Демьянов В. В., Савельева Е. А. и др. Элементарное введение в геостатистику. — М., 1999. — (Пробемы окружающей среды и природных ресурсов / ВИНИТИ; № 11).

Aarts E., Korst J. Simulated Annealing and Boltzmann Machines. — New York: John Wiley & Sons, 1989.

Caers J. Petroleum Geostatistics / Society of Petroleum Engineers. — Richardson, TX, 2005.

Chiles J.-P., Delfiner P. Geostatistics: Modeling Spatial Uncertainty. — New York: John Wiley & Sons, 1999.

Christakos G. Random Field Models in Earth Sciences. — San Diego, CA:

Academic Press, 1992.

Deutsch C. Geostatistical Reservoir Modelling. — [S. l.]: Oxford Univ. Press, 2002.

Deutsch C., Journel A. G. GSLIB: Geostatistical Software Library and User's Guide. — [S. l.]: Oxford Univ. Press, 1998.

Emerly X. Variogram of order : A tool to validate a bivariate distribution model // Mathematical Geology. — 2005. — Vol. 37, N 2. — Р. 163—181.

Goovaerts P. Geostatistics for Natural Resources Evaluation. — [S. l.]:

Oxford Univ. Press, 1997. — 483 p.

Isaaks E. H., Srivastava R. M. An Introduction to Applied Geostatistics. — Oxford, Oxford Univ. Press, 1989.

Journel A. G. Modeling uncertainty: some conceptual thoughts // Geostatistics for the Next Century / Ed. R. Dimitrakopoulos. — Dordrecht:

Kluwer Academic Pub., 1994. — P. 30—43.

Journel A. G., Huijbregts C. J. Mining Geostatistics. — London: Academic Press, 1978. — 600 p.

Mantoglou A., Wilson J. Simulation of random fields with the turning band method / Department of Civil Engineering, M.I.T. — [S. l.], 1981. — (Technical Report N 264).

Metropolis N., Rosenbluth A., Teller A., Teller E. Equations of state calculations by fast computing machines // J. of Chem. Physics. — 1985. — Vol. 21, N 6. — Р. 1087—1092.

Стохастическое моделирование пространственной неопределенности Perrin O., Iovleff S. Estimation a non-stationary spatial structure using simulated annealing // GeoComputation 99 // http://www.geovista.psu.edu/ geocomp/geocomp99/Gc99/028/gc_028.htm.

Soares A. Direct Sequential Simulation and Cosimulation // Mathematical Geology. — 2001. — Vol. 33, N 8. — Р. 911—926.

Tran T. Improving variogram reproduction on dense simulation grids // Computers and Geosiences. — 1994. — Vol. 20. — Р. 1161—1168.

Глава Последовательный геостатистический анализ данных:

примеры исследования 9.1. Использование обычного кригинга для мониторинга радиационного загрязнения в режиме реального времени В данном примере описаны результаты участия обычного кригинга в международном конкурсе сравнения методов пространственной интерполяции (Spatial Interpolation Comparison — SIC2004), организованном Исследовательским центром в Испре, Италия (Joint Research Centre, Ispra, Italy). Подробное описание данных, условий и результатов конкурса опубликовано в специальном выпуске журнала «Applied GIS» в 2005 г. Одним из условий конкурса было использование метода в автоматическом режиме, т. е. настраиваемые параметры метода должны были быть высланы организаторам за сутки до предоставления исходных данных. Для настройки параметров можно было использовать данные по измерению той же величины на той же сети мониторинга, но в другое время.

В конкурсе использовались данные радиационного мониторинга воздуха в районе действующей АЭС, расположенной в Европе [Dubois, Galmarini, 2005].

Таким образом, можно считать, что этот пример демонстрирует возможность применения обычного кригинга для анализа данных мониторинга радиационной обстановки вокруг радиационно опасного объекта в режиме реального времени (on-line) [Savelieva, 2005].

Реальный мониторинг представлен 1008 датчиками. Для проверки качества работы методов в рамках конкурса были выделены тренировочный и валидационный наборы. Тренировочный набор представлен 200 точками, 808 использовались для валидации. Пространственное распределение тренировочных и валидационных точек приведено на рис. 9.1. Сам объект находился в точке с координатами (0, 0).

Последовательный геостатистический анализ данных: примеры исследования Рис. 9.1. Пространственное распределение тренировочных ( ) Использование в режиме реального времени предполагает полностью автоматическое функционирование, т. е. такие параметры метода, как модель пространственной корреляционной структуры и область поиска (область оценки), считаются заданными априори. В данном случае для настройки параметров использовалась историческая информация, 10 наборов измерений того же параметра в 200 тренировочных точках.

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

Разница между максимальным и минимальным значениями в каждой точке не превышала 40, больше 30 разница была только в 6 точках. Визуальную схожесть можно наблюдать на рис. 9.2, где визуализированы данные двух наборов (рис. 9.2а,б) и среднего по 10 наборам (рис. 9.2в).

В. В. Демьянов, Е. А. Савельева Геостатистика: теория и практика Рис. 9.2. Исторические данные для настройки параметров:

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

• Экспериментальные вариограммы для различных наборов измерений также схожи (пример четырех вариограммных изолиний представлен на рис. 9.3): они не обладают анизотропией до расстояния 60 км и демонстрируют анизотропию (больший радиус корреляции в направлении запад—восток) на расстояниях до 200 км. Значение плато и размер области корреляции одинаковы.

• Вариограмма, усредненная по 10 экспериментальным вариограммам, отражает все свойства отдельных вариограмм (рис. 9.4а). Можно предположить, что она будет отражать пространственную корреляционную структуру и других наборов измерений этой переменной.

По усредненной вариограмме построена модель, которая и используется в рамках обычного кригинга (рис. 9.4б). Выбрана сферическая модель со значением наггета 33,59, плато — 267,0 и эллипсом корреляции с радиусами 306,3 и 230,4 м, главная ось — в направлении запад—восток.

Зона поиска представлена эллипсом с радиусами 310 и 235 м, главная ось — в направлении запад—восток.

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

Для картирования были предоставлены два набора данных: обычный выброс (данные, статистически аналогичные априорным данным) и аварийный выброс (искусственно смоделированный выброс, наложенный на обычные Последовательный геостатистический анализ данных: примеры исследования данные) [Dubois, Galmarini, 2005]. Обычный кригинг с описанными выше параметрами был использован для обоих наборов данных при прогнозировании значений в 808 валидационных точках. Полученные обычным кригингом оценки представлены на рис. 9.5. Результаты интерполяции различаются, так как зависят не только от параметров модели, но и от исходных данных.

Рис. 9.3. Примеры экспериментальных вариограмм Рис. 9.4. Усредненная по 10 экспериментальным вариограммам В. В. Демьянов, Е. А. Савельева Геостатистика: теория и практика Рис. 9.5. Результат интерполяции обычным кригингом для обычных данных (а) Как и ожидалось, обычный набор статистически похож на исторические наборы, поэтому результат валидации соответствует результатам кроссвалидации (коэффициент корреляции — 0,78).

Рис. 9.6. Моделирование обычным кригингом экстремальных значений Пятно, характеризующее выброс, обнаружено и хорошо видно на рисунке.

Но другие области (особенно северная часть) кажутся не подверженными влиянию выброса и выглядят одинаково для обоих исходных наборов данных. Пятно растянуто и сглажено (рис. 9.6), но обнаружено. Основная проблема состоит в том, что в окрестности вокруг выброса (зоне максимальноГлава Последовательный геостатистический анализ данных: примеры исследования го градиента) нарушено предположение обычного кригинга о постоянстве среднего. Это вызывает искажение оценки (светлые пятна вокруг выброса на рис. 9.5б). Таким образом, оценка кригинга вокруг выброса в момент выброса не может считаться корректной.

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

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

Рис. 9.7. Вариация кригинга одинаковая для обычных данных Воспроизведение статистических характеристик на валидационном наборе приведено в табл. 9.1. Обычный набор характеризуется хорошим совпадением. Для набора с выбросом хорошо воспроизведены среднее и медиана.

Недооценка максимума и стандартного отклонения соответствует сглаженности оценки. Значение минимума в кригинговой оценке меньше минимума исходных данных, что связано с искажением оценки в области максимального градиента переменной. Там обнаружено 6 таких точек.

В. В. Демьянов, Е. А. Савельева Геостатистика: теория и практика Таблица 9.1. Сравнение статистических характеристик, оцененных обычным Реальные обычные Оцененные обычные Реальные с выбросом Оцененные с выбросом Валидационные ошибки анализировались с помощью таких характеристик, как среднее от абсолютных значений ошибок (САО), среднее ошибок (СО), коэффициент корреляции между оцененными и реальными значениями и корень из среднеквадратичной ошибки (СКО). Эти характеристики представлены в табл. 9.2. Результаты прогноза для обычных данных лучше, чем для данных с выбросом, что вполне очевидно, так как обычные данные по характеристикам совпадают с теми, по которым настраивались параметры модели. Тем не менее результаты по данным с выбросом не являются бессмысленными.

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

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

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

Последовательный геостатистический анализ данных: примеры исследования 9.2. Анализ неопределенности в моделировании гидрогеологической структуры Этот пример описывает моделирование одного гидрогеологического осадочного слоя в рамках гидрогеологической системы из 10 слоев. Анализ данных проводился в рамках совместных исследований ИБРАЭ РАН и Pacific Northwest National Laboratory по программе РАН и Министерства энергетики США. Результаты исследований представлены в [Savelieva et al., 2002].

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

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

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

Моделирование проводилось последовательно для каждого слоя. Здесь рассмотрен один из слабо проводящих слоев (U4), являющийся очень важным в данной гидрогеологической системе.

Моделирование проводилось в два этапа:

• получение зоны присутствия данного слоя (задача бинарной классификации);

• оценка толщины гидрогеологического слоя в областях его присутствия.

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

Исходный набор данных содержал 401 скважину, где измерялись толщины гидрогеологических слоев. На рис. 9.8 приведено пространственное расположение скважин, каждая помечена в соответствии с присутствием (отсутствием) в ней слоя U4. Несколько скважин обозначены как неопределенные, т. е. эксперт-геолог не смог окончательно решить, присутствует В. В. Демьянов, Е. А. Савельева Геостатистика: теория и практика ли в керне слой U4. Толщина слоя U4 измерена в 117 скважинах, где он определенно обнаружен.

Рис. 9.8. Пространственное распределение присутствия Первая часть задачи — бинарная классификация — может решаться с помощью индикаторного кригинга. Индикаторное преобразование:

Индикаторный кригинг даст вероятность присутствия слоя U4.

Но сначала требуется провести построение и моделирование индикаторной вариограммы. На рис. 9.9 приведены экспериментальная индикаторная вариограмма и ее модель. Они визуализированы с помощью изолиний вариограммы. Параметры модели вариограммы: наггет — 0,1, плато — 0,16, радиусы корреляции — 29 030 и 10 480 м, главная ось эллипса — в направлении 30° по часовой стрелке от оси север—юг.

Последовательный геостатистический анализ данных: примеры исследования Рис. 9.9. Экспериментальная индикаторная вариограмма присутствия слоя U4 (а) Результат кросс-валидации дает ошибку классификации 18%.

Расчет проводился на ячейках размером 150150 м и на области, ограниченной рассматриваемой гидрогеологической моделью. Результат картирования вероятности присутствия слоя U4 представлен на рис. 9.10.

Рис. 9.10. Вероятность присутствия гидрогеологического слоя В. В. Демьянов, Е. А. Савельева Геостатистика: теория и практика Второй шаг — стохастическое моделирование толщины гидрогеологического слоя U4 с учетом его присутствия. Присутствие определяется по вероятности, полученной на предыдущем шаге. Эту вероятность можно учитывать по-разному. Можно сразу провести классификацию, считая, что слой присутствует, если вероятность его присутствия больше пороговой вероятности (например, 0,5). Другой подход состоит в розыгрыше присутствия каждый раз при движении по сетке оценивания. Мы использовали оба варианта.

В качестве метода стохастического моделирования использовалось последовательное гауссово моделирование. Данные подвергались нормализующему преобразованию. Бинормальность проверялась эмпирическим тестом — с использованием мадограммы:

Результат теста приведен на рис. 9.11.

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

Примеры полученных реализаций приведены на рис. 9.12 (на рис. 9.12б — без розыгрыша присутствия слоя U4). В табл. 9.3 собраны некоторые статистические характеристики для пяти произвольно выбранных реализаций.

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

Последовательный геостатистический анализ данных: примеры исследования Рис. 9.12. Примеры реализаций толщины гидрогеологического слоя Таблица 9.3. Глобальные статистические характеристики В. В. Демьянов, Е. А. Савельева Геостатистика: теория и практика Рис. 9.13. Экспериментальные вариограммы Аналогично можно проводить моделирование для других гидрогеологических слоев.

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

Воспользуемся данными по радиоактивному загрязнению почвы 241Am, которые использовались в рамках совместных исследований ИБРАЭ РАН и Sandia National Laboratories по программе РАН и Министерства энергетики США. Результаты исследований опубликованы в [Kanevski et al., 2006].

Данные представляют собой набор измерений гамма-детектором в ряде точек, покрывающих большую площадь. Исходные 193 измерения были использованы для моделирования пространственного поля загрязнения и оценки в 917 валидационных точках. Значения в валидационных точках были изначально скрыты от исследователей для обеспечения чистоты эксперимента и приведены лишь после получения оценок для сравнения качества моделирования каждым методом. Исходные измерения приведены на рис. 9.14. Валидационные значения представлены на рис. 9.15.

Целью исследования был вероятностный прогноз превышения уровней загрязнений 17, 27 и 38 пКи/г.

Последовательный геостатистический анализ данных: примеры исследования Рис. 9.14. Исходные данные по загрязнению 241Am в 193 точках Рис. 9.15. Валидационные данные по загрязнению 241Am в 917 точках Одним из ключевых факторов успеха при решении проблемы валидации является репрезентативность исходных данных. Так, если исходные данные сильно отличаются от валидационных, трудно качественно оценить значения в валидационных точках на основе исходных данных. Исходный и валидационный наборы данных достаточно однородно распределены в пространстве. Глобальные статистики для обоих наборов, приведенные в табл. 9.4, имеют близкие значения, из чего следует хорошая репрезентативность исходных данных. Заметим, однако, что длинный хвост высоких значений для валидационных данных вдвое превышает максимальное значение исходных данных.

В. В. Демьянов, Е. А. Савельева Геостатистика: теория и практика Таблица 9.4. Итоговая статистика для данных по загрязнению Am214 (пКи/г) Дальнейшее сравнение исходных и валидационных наборов заключается в рассмотрении их пространственной корреляции. Сравнение вариограмм по всем направлениям для исходных и валидационных данных показывает их близость за исключением более высокой вариабельности и меньшей стационарности валидационных данных (рис. 9.16). При рассмотрении вариограмм по различным направлениям можно выявить более значительные различия между корреляционными структурами исходных и валидационных данных. Так, валидационные данные имеют четкую геометрическую анизотропию в горизонтальном направлении восток—запад, в то время как исходные данные демонстрируют слабую анизотропию в вертикальном направлении север—юг (рис. 9.17а). Вариограммная модель, построенная на основе исходных данных, приведена на рис. 9.17б. Вариограммные модели были также построены для преобразованных значений — нормализованных и индикаторных переменных для 9 пороговых значений (для использования в соответствующих методах). Качество всех моделей было проверено при помощи кросс-валидации и тестирования на 30 данных, предварительно отделенных от исходного набора.

Последовательный геостатистический анализ данных: примеры исследования Рис. 9.16. Вариограммы по всем направлениям для исходных (нижняя линия) Рис. 9.17. Изолинии вариограммной розы для исходных сырых данных (а) Разнообразные геостатистические модели были применены для решения задачи валидации. Использовались геостатистические оцениватели: простой (SK) и обычный (OK) кригинг для получения точечных оценок загрязнения, а также индикаторный кригинг (IK) для получения вероятностных оценок. Три стохастические модели — последовательное гауссово моделирование (SGS), последовательное индикаторное моделирование (SIS) и моделирование отжига (SA). Таким образом, сравнительный анализ методов был проведен для широкого спектра геостатистических моделей. Расчеты простого и обычного кригинга проводились при помощи пакета программ «Геостат Офис» [Kanevski, Maignan, 2004], результаты стохастического моделирования и индикаторного кригинга были получены при помощи программ GSLIB [Deutsch, Journel, 1998].

В. В. Демьянов, Е. А. Савельева Геостатистика: теория и практика Таблица 9.5. Сравнение итоговой статистики для валидационных данных Нижний квартиль (25%) Верхний квартиль (75%) Стандартное отклонение Примечание. SK — простой кригинг, OK — обычный кригинг, IK — усредненная оценка индикаторного кригинга, SGS — среднее значение гауссова моделирования, SIS — среднее значение индикаторного моделирования, SA — среднее значение моделирования отжигом. Полужирным шрифтом выделены значения статистик оценок, наиболее близкие к статистикам валидационного набора.

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

Индикаторным кригингом с использованием девяти индикаторных переменных были получены локальные функции распределения вероятности в валидационных точках. Для сравнения с оценками кригинга были использованы усредненные значения Е-типа (E-type) (см. Раздел 7.2).

Стохастическое моделирование было проведено на регулярной сетке с шагом 55 м. Далее значения в валидационных точках были получены методом ближайшего соседа. Такая последовательность обусловлена ограничениями программ пакета GSLIB, что может внести смещение в окончательные результаты валидации. Однако выбор достаточно высокого разрешения сетки моделирования по сравнению с разрешением сети валидационных данных позволяет считать это хорошей аппроксимацией. Ошибка аппроксимации при использовании в данном случае метода ближайшего соседа Последовательный геостатистический анализ данных: примеры исследования будет значительно меньше, чем ошибка измерения и локальная неопределенность в валидационных точках. При анализе результатов стохастического моделирования вместо использования усредненных оценок Е-типа (см. Главу 8) были рассчитаны статистические параметры распределений отдельных реализаций. Далее эти статистики были усреднены для сравнения с валидационным распределением. Как видно из табл. 9.5, усредненная оценка индикаторного кригинга (IK) сильно сглажена по сравнению с валидационными данными. Это обусловлено выбором среднего значения промоделированных распределений. Оценка кригинга дала хорошее совпадение медианы распределения и, более того, эксцесса, который близок к экцессу валидационного распределения. Стохастическое моделирование позволяет лучше воспроизвести валидационное распределение, чем кригинг. Можно видеть, что различные статистические параметры воспроизводятся лучше разными методами. Так, результаты гауссова моделирования (SGS) имеют значения минимума, максимума и медианы, наиболее близкие к соответствующим параметрам валидационного набора. Индикаторное моделирование (SIS) дало наиболее близкие средние значения и значения квартилей. Реализации моделирования отжигом (SA) имеет наилучшие стандартное отклонение и коэффициент симметрии. Еще раз подчеркнем, что статистические параметры были получены путем усреднения статистик каждой из 100 реализаций для каждого алгоритма.

Пространственные корреляционные структуры оценок кригинга и реализаций стохастического моделирования представлены соответствующими вариограммами, которые сравнивались с вариограммой валидационных данных. Вариограмма оценки кригинга ожидаемо недооценивает уровень пространственной вариации, хотя размеры корреляции представлены достаточно хорошо (рис. 9.18д). Вариограмма оценки кригинга не имеет наггета, в то время как наггет по вариограмме валидационных данных составляет 25—30% априорной вариации.

Стохастическое моделирование позволяет лучше промоделировать вариабельность и неопределенность пространственной корреляции, которую можно представить доверительными интервалами вариограммы на основе стохастических реализаций. Распределение вариограмм для реализаций представлено средней вариограммой с доверительным интервалом ±2 (рис. 9.19). Как было отмечено выше, пространственная корреляция валидационных данных существенно отличается от корреляции исходных данных. Все алгоритмы кроме индикаторного моделирования (SIS) демонстрируют хорошее совпадение пространственной корреляционной структуВ. В. Демьянов, Е. А. Савельева Геостатистика: теория и практика ры. Усредненная вариограмма реализаций индикаторного моделирования (SIS) недооценивает уровень вариабельности даже с учетом интервала неопределенности (см. рис. 9.19б). Широкие доверительные интервалы вариограмм реализаций говорят о значительной неопределенности, которую воспроизводят стохастические реализации. Усредненные вариограммные розы реализаций демонстрируют хорошее совпадение со структурой валидационных данных (см. рис. 9.18). Результаты моделирования отжига (SA) дают наиболее близкое совпадение пространственной корреляции с валидационным распределением (см. рис. 9.18а). Усредненные вариограммные розы для реализаций гауссова (SGS) и индикаторного (SIS) моделирования более близки к пространственной корреляции исходного распределения (см. рис. 9.18б,г).

Рис. 9.18. Сравнение контуров вариограммной розы валидационных данных с усредненными вариограммами по 100 стохастическим реализациям гауссова моделирования невязок нейронной сети (г), Последовательный геостатистический анализ данных: примеры исследования Рис. 9.19. Средняя вариограмма по всем направлениям с доверительным интервалом ±2 для 100 реализаций:

гауссово моделирование невязок нейронной сети (г) в сравнении с вариограммой для валидационных данных Анализ невязок, оставшихся после оценок кригинга, обнаружил существование пространственной корреляции в них на малых расстояниях. Возможно, это связано с присутствием нестационарности, которая не учитывается кригингом. Так, крупномасштабный пространственный тренд, выявленный в валидационных данных, не проявлялся в исходных данных и поэтому не был учтен. Оценки обычного кригинга (OK) вместе c соответствующей вариацией приведены на рис. 9.20.

Рис. 9.20. Оценки обычного кригинга (а) и соответствующая вариация оценки (б) В. В. Демьянов, Е. А. Савельева Геостатистика: теория и практика Вероятностные оценки превышения уровня загрязнения 241Am 27 пКи/г, полученные индикаторным кригингом (IK) и различными алгоритмами стохастического моделирования, приведены на рис. 9.21.

Рис. 9.21. Оценки вероятности превышения уровня концентрации Am241 27 пКи/г, полученные гауссовым моделированием (SGS) (а), индикаторным моделированием (SIS) (б), моделированием отжига (SA) (в), индикаторным кригингом (IK) (г) Качество валидационной оценки можно определить путем сравнения промоделированных локальных функций распределения с валидационными данными. На рис. 9.22 представлены локальные функции распределения, полученные различными методами в четырех валидационных точках. На те же графики нанесены оценки обычного кригинга (OK) с доверительным интервалом ± ошибки оценки кригинга.

Количественный анализ удовлетворения доверительных интервалов, построенных на основе локальных функций распределения во всех 917 валидационных точках, приведен в табл. 9.6. Были выбраны четыре доверительных интервала: размах (между минимумом и максимумом оценки), 90% распределения между 5 и 95%; 80% распределения между 10 и 90% распределения; межквартильный интервал между 25 и 75% распределения. Из табл. 9.6 видно, что все три метода стохастического моделирования дают доверительные интервалы, хорошо удовлетГлава Последовательный геостатистический анализ данных: примеры исследования воряющие валидационным данным. Процент валидационных данных, попавших в соответствующий интервал, хорошо согласуется с размерами интервала. Наилучшие результаты показало моделирование отжига (SA).

cdf cdf Рис. 9.22. Локальные кумулятивные функции распределения вероятности в четырех валидационных точках, полученные различными методами: индикаторным кригингом (IK), гауссовым моделированием (SGS), индикаторным моделированием (SIS), моделированием отжига (SA), моделированием невязок нейронной сети (MLRSGS), в сравнении с соответствующими валидационными измерениями (жирная вертикальная линия) и значениями оценки обычного кригинга (OK — тонкая вертикальная линия) c доверительным интервалом ± Таблица 9.6. Доля валидационных данных, попавших в интервалы неопределенности, полученные различными методами стохастического моделирования, % В. В. Демьянов, Е. А. Савельева Геостатистика: теория и практика Пример исследования показал, что методы стохастического моделирования способны хорошо воспроизводить пространственную вариабельность.

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

Литература Deutsch C. V., Journel A. G. GSLIB: Geostatistical Software Library and User’s Guide. — New York: Oxford Univ. Press, 1998. — 369 p.

Dubois G., Galmarini S. Introduction to the Spatial Interpolation Comparison (SIC) 2004 exercise and presentation of the data sets // Applied GIS. — 2005. — Vol. 1, N 2. — Р. 9.1—9.10 (http://publications.epress.monash.edu/ doi/pdf/10.2104/ag050009).

Kanevski M., Demyanov V., Savelieva E. et al. Validation Of Geostatistical And Machine Learning Models For Spatial Decision-Oriented Mapping // Proceeding of StatGIS 99 / Ed. J. Piltz, J. Heyn. — Klagenfurt, 2006.

Kanevski M., Maignan M. Analysis and modelling of spatial environmental data. — Lausanne: EPFL Press, 2004. — 288 p. — (With a CD and educational/research MS Windows software tools).

Savelieva E. Using Ordinary Kriging to Model Radioactive Contamination Data // Applied GIS. — 2005. — Vol. 1, N 2. — Р. 10.1—10.10.

Savelieva E., Kanevski M., Timonin V. et al. Uncertainty in the hydrogeologic structure modeling // Proceedings of IAMG2002 conference. — [S. l.], 2002. — Р. 481—486.

Глава Комбинированные модели ИНС и геостатистики Раздел 10.1 настоящей главы посвящен постановке задачи комбинированного моделирования на основе искусственных нейронных сетей (ИНС) и геостатистики. В Разделах 10.2, 10.3 приведены примеры использования предложенного комбинированного метода для моделирования пространственных и временных данных: в Разделе 10.2 рассмотрено картирование атмосферных осадков при помощи кригинга невязок ИНС [Kanevsky et al., 1998], а в Разделе 10.3 — прогнозирование электропотребления при помощи стохастического моделирования невязок ИНС.

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

10.1. Геостатистический анализ невязок Идея метода заключается в моделировании нелинейного крупномасштабного тренда при помощи ИНС и последующего моделирования невязок геостатистическими методами. Этот подход был впервые предложил М. Ф. Каневский В. В. Демьянов, Е. А. Савельева Геостатистика: теория и практика [Kanevsky et al., 1996а] и в дальнейшем развивался в последующих работах [Kanevsky et al., 1997, 1998; Demyanov et al., 2000; 2001; Savelieva et al., 2000]. Изначально невязки, оставшиеся после применения ИНС для пространственного оценивания данных, были оценены обычным кригингом.

Таким образом, итоговая оценка была получена как сумма оценки ИНС и кригинга невязок. Последующее развитие подхода касалось обобщения на случай нескольких переменных и применения стохастического моделирования невязок. В случае нескольких переменных используются ИНС с несколькими выходными нейронами и для моделирования невязок применяется кокригинг (см. Главу 6) [Kanevski et al., 1997]. Методы стохастического моделирования (см. Главу 8) привлекаются для оценивания невязок аналогичным образом [Demyanov et al., 2000]. В Разделе 9.3 был приведен один из таких методов с использованием последовательного гауссового моделирования в качестве сравнения (MLRSGS). Комбинированный подход на основе ИНС и геостатистики нашел применение и в других работах, например в [Cortez et al., 1998; Bryan, Adams, 2002; Zhang et al., 2004].

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

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

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

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

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

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

Кригинг невязок (residual kriging), как и универсальный кригинг, предполагает, что неизвестное среднее значение m(x) меняется во всей области исследования S так, что нельзя допустить постоянство даже локального среднего. В этом случае компонента тренда моделируется отдельно, используя другие математические или физические подходы. В качестве модели тренда можно использовать прогноз, выполненный с использованием физической модели процесса, являвшегося причиной формирования поля Z(x). Компонента тренда, как уже указывалось, может моделироваться с помощью нелинейных адаптивных методов (например, искусственных нейронных сетей, вейвлетов [Demyanov et al., 2001], метода регрессии на опорных векторах [Kanevski et al., 2003, 2004] и др.), использующих набор измерений как информацию для настройки своих параметров. После выделения тренда простой или обычный кригинг используется на невязках модели к измеренным значениям поля.

В. В. Демьянов, Е. А. Савельева Геостатистика: теория и практика Кригинг невязок можно также интерпретировать как гибридную модель, объединяющую различные по математической или физической базе методы.

10.2. Пример использования кригинга невязок В данном примере использованы метеорологические данные по усредненным за 10 дней выпадениям осадков в Швейцарии в 1986 г. Задача состояла в оценке значений в 367 точках по 100 измерениям (рис. 10.1). Эти данные распространялись в рамках международного конкурса сравнения методов пространственной интерполяции (Spatial Interpolation Comparison — SIC’97), организованного геостатистическим порталом AI-GEOSTAT [SIC’97] и Группой мониторинга радиоактивности в окружающей среде Института окружающей среды в Объединенном исследовательском центре (Radioactivity Environmental Monitoring group of the Environment Institute at the Joint Research Centre), Испра, Италия.

Рис. 10.1. Исходные данные (100 точек), использованные для моделирования (а), и валидационные данные для проверки качества оценки (367 точек) (б) Данные обладают корреляционной структурой с трендом — дрейфом (определение дрейфа см. в Разделе 4.2), который почти во всех направлениях практически монотонно убывает (рис. 10.2). Очевидно, что для картирования таких данных требуется моделирование тренда. Было предложено построить нелинейную модель тренда с помощью искусственной нейронной сети [Kanevsky et al., 1998] или с использованием вейвлет-метода [Demyanov et al., 2001].