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

Упражнение 5.6. Сравнение оценок обычного кригинга и метода обратных квадратов расстояния На рис. 5.4. приведены три интерполяционные оценки на основе пяти данных. Две из них получены обычным кригингом с большим и маленьким радиусами корреляции, третья — методом обратных квадратов расстояния (см. Раздел 3.1). Поставить в соответствие оценкам A, B, C приведенные выше интерполяционные модели.

Рис. 5.4. Интерполяционные оценки A, B, C на основе пяти данных:

обычный кригинг с большим и маленьким радиусами корреляции Пример результата обычного кригинга Для более полного понимания функционирования обычного кригинга рассмотрим пример его применения на искусственных данных, моделирующих пористость. Пространственное распределение исходных данных приведено на рис. 5.5 — шесть точек, имеющих значения в диапазоне от 0,05 до 0,30. Цель примера — показать влияние на результат обычного кригинга модели пространственной корреляционной структуры (вариоГлава Геостатистические интерполяции для одной переменной граммы). При этом будем рассматривать оценку обычного кригинга и вариацию обычного кригинга.

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

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

В. В. Демьянов, Е. А. Савельева Геостатистика: теория и практика Рис. 5.6. Оценка обычного кригинга при изотропной модели пространственной Геостатистические интерполяции для одной переменной Рис. 5.7. Вариация обычного кригинга при изотропной модели пространственной В. В. Демьянов, Е. А. Савельева Геостатистика: теория и практика На рис. 5.8—5.13 представлены оценка и вариация кригинга для случаев с моделью вариограммы с геометрической анизотропией: рис. 5.8 и 5. соответствуют направлению главной оси вдоль оси X, рис. 5.10 и 5.11 — вдоль оси Y, рис. 5.12 и 5.13 — направлениям 45°, 30° и 60° от оси X против часовой стрелки. Радиус корреляции вдоль главной оси равен 40 м. Значение радиуса корреляции вдоль малой оси варьируется. Везде заметно влияние направления главной оси, а изменение значения второго радиуса корреляции приводит к результатам, аналогичным описанным выше для изотропного случая.

Рис. 5.8. Оценка обычного кригинга при анизотропной модели пространственной корреляционной структуры (вариограмме), направление главной оси — вдоль оси X, радиус корреляции вдоль главной оси — 40 м, значения радиуса корреляции вдоль малой оси: а — 5 м; б — 10 м; в — 20 м; г — 30 м Геостатистические интерполяции для одной переменной Рис. 5.9. Вариация обычного кригинга при анизотропной модели пространственной корреляционной структуры (вариограмме), направление главной оси — вдоль оси X, радиус корреляции вдоль главной оси — 40 м, значения радиуса корреляции вдоль малой оси: а — 5 м; б — 10 м; в — 20 м; г — 30 м Рис. 5.10. Оценка обычного кригинга при анизотропной модели пространственной корреляционной структуры (вариограмме), направление главной оси — вдоль оси Y, радиус корреляции вдоль главной оси — 40 м, значения радиуса корреляции вдоль малой оси: а — 10 м; б — 20 м; в — 30 м В. В. Демьянов, Е. А. Савельева Геостатистика: теория и практика Рис. 5.11. Вариация обычного кригинга при анизотропной модели пространственной корреляционной структуры (вариограмме), направление главной оси — вдоль оси Y, радиус корреляции вдоль главной оси — 40 м, значения радиуса корреляции Рис. 5.12. Оценка обычного кригинга при анизотропной модели пространственной корреляционной структуры (вариограмме), Геостатистические интерполяции для одной переменной Рис. 5.13. Вариация обычного кригинга при анизотропной модели пространственной корреляционной структуры (вариограмме), радиус корреляции вдоль главной оси — 40 м, радиус корреляции вдоль малой оси — 20 м, направление главной оси: а — 45°, б — 30°; в — 60° (от оси X против часовой стрелки) На рис. 5.14 и 5.15 показано влияние на оценку и вариацию кригинга вариации наггета. Во всех рассмотренных до этого случаях значение данного параметра было равно нулю. Здесь рассматривается случай с изотропной моделью пространственной корреляции, радиус корреляции равен 5. При таком радиусе корреляции большая часть оцениваемой области соответствует среднему значению оценки. Это усредненное значение больше всего подвержено влиянию наггета. Рассмотрены случаи со значением наггета 0, и 0,9. Картина оценки при этом меняется не очень существенно, но можно заметить, что зоны низких и высоких значений (вокруг точек) уменьшаются, как при уменьшении радиуса корреляции. А значение вариации кригинга существенно возрастает при увеличении этого параметра (можно сравнить шкалы на рис. 5.15а и 5.15б). Наггет характеризует уровень неопределенности модели пространственной корреляции.

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

В. В. Демьянов, Е. А. Савельева Геостатистика: теория и практика Рис. 5.14. Оценка обычного кригинга при изотропной модели пространственной корреляционной структуры (вариограмме), радиус корреляции 5 м, Рис. 5.15. Вариация обычного кригинга при изотропной модели пространственной корреляционной структуры (вариограмме), радиус корреляции 5 м, 5.4. Универсальный кригинг Универсальный кригинг, или кригинг с трендом (universal kriging, UK), предполагает, что неизвестное среднее значение m(x) плавно меняется во всей области исследования S. В некоторых случаях невозможно предположить локальное постоянство среднего даже в окрестности оцениваемой точки W(x). Одним из возможных в таком случае подходов является именно универсальный кригинг. Предполагается, что детерминистическая компонента случайной переменной (тренд) моделируется как линейная комбинация K + 1 базисных (известных) функций fk(x) (по принятому соглашению f 0 ( x ) = 1) с коэффициентами ak(x), неизвестными и постоянными внутри окрестности оцениваемой точки x W(x):

Рассмотрим, как выполнить в таком случае условие несмещенности оценки:

Получаем набор из K + 1 дополнительных ограничений, но избавляемся от необходимости оценивать коэффициенты ak(x):

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

В. В. Демьянов, Е. А. Савельева Геостатистика: теория и практика Вариация универсального кригинга также может быть записана следующим Здесь следует обратить внимание, что в системе уравнений универсального кригинга используются ковариации CR() для случайной компоненты R(x) функции Z(x), априорное знание которых предполагается. Кроме того, требуется априорное знание набора базисных функций. Подробнее ознакомиться с тем, как на практике подходят к решению задач подготовки модели тренда и ковариационной функции остатков, можно, например, в [Goovaerts, 1997; Armstrong, 1984]. Но, вообще говоря, универсальный кригинг не получил широкого распространения, так как задача подбора функций для моделирования тренда не является прозрачной.

5.5. Логнормальный кригинг Логнормальным случайным процессом {Z ( x ) : x S } называется такой положительно-определенный процесс, когда Y(x) является гауссовым процессом:

Пусть Y*(x0) — оценка функции Y(x) в точке x0, где нет измерения, полученная с помощью кригинга на основании известных данных Y ( x1 ), Y ( x2 ),..., Y ( xn ) в точках измерений x1, x2,..., xn. Теперь требуется получить оценку функции Z(x) в этой точке. Если получать такую оценку, просто делая обратное логарифму преобразование то она будет смещенной [Dowd, 1982], т. е.

А цель любого кригинга, в том числе и логнормального, — наилучшая несмещенная оценка. Поэтому, чтобы учесть и исправить возникающее при обратном преобразовании смещение, используется формула Здесь вариация Var Y ( x0 ) представляет собой не вариацию обычного кригинга OK ( x0 ), а значение, определяемое как Y = [Y ( x1 ), Y ( x2 ),..., Y ( xn )]. Таким образом, несмещенная оценка Z(x0) может быть получена на основании следующей формулы:

где — множитель Лагранжа, значение которого находится при решении системы уравнений обычного кригинга для оценки Y*(x0).

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

Пример данных для логнормального кригинга Примером могут служить данные по пространственному распределению крабов. Они получены от Всероссийского научно-исследовательского института рыбного хозяйства и океанографии (ВНИРО) и представляют собой результат траловой съемки краба Берди в Беринговом море в 2003 г. Для крабов характерны огромные скопления на фоне практически нулевых значений. Разброс значений на пять порядков не дает возможности оценивать пространственную корреляцию с помощью вариограммы (рис. 5.16).

В. В. Демьянов, Е. А. Савельева Геостатистика: теория и практика Рис. 5.16. Экспериментальная вариограмма для данных траловой съемки пространственного распределения краба Берди После логарифмического преобразования данные вполне пригодны для геостатистического анализа и моделирования — вариограмма на рис. 5. легко моделируется сферической моделью.

Рис. 5.17. Экспериментальная вариограмма для данных траловой съемки пространственного распределения краба Берди Логарифмически преобразованные данные не всегда соответствуют требованию нормальности. Но невыполнение этого условия нарушает корректность обратного преобразования (5.22), что может привести к сомнительной по качеству оценке. В случае с крабами корректность соблюдается, их распределение отлично соответствует логнормальному (рис. 5.18).

Рис. 5.18. График соответствия эмпирического распределения теоретическому (probability-probability plot) для данных траловой съемки пространственного Геостатистические интерполяции для одной переменной 5.6. Некоторые дополнительные аспекты кригинга 5.6.1. Веса кригинга и эффект экранирования Выше (при обсуждении основных свойств простого кригинга — см. раздел 5.2) упоминалось, что веса кригинга не зависят от известных значений функции, а определяются моделью пространственной корреляции. Отметим еще несколько свойств весов кригинга, которые позволят лучше понять кригинг как метод.

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

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

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

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

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

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

Все описанные выше свойства можно проверить. Для этого достаточно рассмотреть точки 1—5 на рис. 5.19 как пространственную модель точек для оценки в центре (большой серой точке, помеченной «?»). В качестве ковариационной функции можно, например, использовать изотропную сферическую модель с нулевым наггетом, единичным плато и радиусом корреляции 1 км.

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

Геостатистические интерполяции для одной переменной 5.6.2. Вариация кригинга и неопределенность оценки Модели кригинга позволяют получить оценки локального среднего (оценка кригинга) и локальной вариации (вариация кригинга). Полученную вариацию кригинга можно использовать для описания неопределенности оценки. Если принять гипотезу о мультинормальности случайных переменных, то 95%-ные доверительные интервалы, в которых лежит истинное значение функции Z(x0), будут определяться как где — стандартное отклонение, полученное из кригинговой вариации.

Если не делать предположения о мультинормальности, то размер 95%-ного доверительного интервала увеличивается [Chiles, Delfiner, 1999] (отсутствие знания увеличивает неопределенность) до 6. Это прямое следствие неравенства Высочанского — Петунина, полученного в 1980 г. Единственное предположение является достаточно слабым, предполагается, что распределение ошибки является непрерывным и унимодальным.

Неравенство Высочанского — Петунина формулируется следующим образом: если X — случайная переменная с плотностью распределения f, неубывающей до моды v, а потом невозрастающей, и если — ожидаемое стандартное отклонение от произвольного значения, то При сформулированных предположениях доверительный интервал является примерно 90%-ным. Чтобы получить 95%-ный интервал, требуется положить t = 3, так как 4 = 0, 049.

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

Изолиния оценки будет с заданной вероятностью (например, 0,9) лежать между контурами доверительных интервалов, которые образуют «кориВ. В. Демьянов, Е. А. Савельева Геостатистика: теория и практика дор». Ширина «коридора» характеризует неопределенность изолинии оценки. Эта так называемая «толстая» изолиния содержит в себе с определенной вероятностью изолинию действительного распределения, истинное положение которой в точности неизвестно. Пример «толстой» изолинии, представляющей 90%-ный доверительный интервал (±2ОК), приведен на рис. 5.20. Этот пример относится к анализу данных по загрязнению поверхности в результате Чернобыльского выброса в мае 1986 г. [Kanevski et al., 1999]. Результаты измерений картированы с использованием обычного кригинга. Изолиния относится к одному из критических уровней загрязнения — 1000 Ки/км2.

Рис. 5.20. «Толстые» изолинии, представляющие 90%-ный доверительный интервал для изолинии 1000 Ки/км2 по результатам картирования обычным кригингом загрязнения поверхности 137Cs в районе Чернобыльского выброса С другой стороны, в соответствии с (5.14) для простого кригинга и (5.18) или (5.19) для обычного кригинга, вариация кригинга не зависит от значений исходных данных и зависит от ковариационной функции и конфигурации (взаимного расположения) данных. Это также было отмечено и в искусственном примере на обычный кригинг. Таким образом, кригинг дает одинаковую вариацию для всех точек с аналогичной конфигурацией исходных данных (при использовании одной и той же модели пространственной корреляционной структуры) и для случаев, когда исходные данные близки по значениям и когда они сильно различаются, т.е. фактически вариация Геостатистические интерполяции для одной переменной кригинга является характеристикой плотности исходных данных. Это можно наблюдать на рис. 5.21, где изображена вариация кригинга вместе с исходными данными, полученная для оценки кригинга. Видно, что вариация оценки не зависит от значения самой оценки. Зависимость вариации кригинга только от плотности исходных данных позволяет использовать ее для оптимизации сети мониторинга.

(крестиками отмечены точки расположения исходных данных) 5.6.3. Учет собственных ошибок измерений в уравнениях кригинга Бывают случаи, когда измеренные значения сопровождаются ошибкой измерения, полученной, например, как сведения о систематической ошибке прибора или методологии измерений, т.е. вместе со значением функции Z(xi) имеется еще и ошибка измерения, представленная как локальная вариация (xi) = i. Если сами ошибки не являются пространственно коррелированными, то система уравнений обычного кригинга может быть видоизменена так, чтобы учитывать ошибки при оценке [Kanevski et al., 1993]:

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

Увеличение неопределенности (кригинговой вариации) вызвано введением дополнительной неопределенности в данные. Обычный кригинг предполагает все данные абсолютно точными.

5.6.4. Блочный и точечный кригинг Любое измерение Z(x) всегда соотносится с некоторым ненулевым конечным объемом. Это может быть кусок породы или почвы, где берется проба.

Обычно размер измерения позволяет приписать его к точке с координатой x. Размер оцениваемого значения такой же, как и измеряемого, и оценка также приписывается к координате. Но иногда целью оценки является среднее по определенному объему, если именно такой размер соответствует определенным действиям (например, очистке и т. п.). Блочный кригинг — обобщенное название метода для определения среднего значения функции z по какому-либо измеримому сегменту (длине, площади, объему) любого размера или формы в противоположность точечному кригингу, относящемуся к нулевому размеру оценки (пробы).

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

Система уравнений блочного кригинга выглядит точно так, как система уравнений обычного кригинга (5.17). Единственное различие состоит в том, что в правой части системы стоит ковариация не для двух точек, а для блока (объема) и точки C ( x,V ( x ) ). Определяется такая ковариация в соответствии с формулой [Journel, Huijbregts, 1978] На практике такая ковариация определяется как среднее точечных ковариаций точки x и N точек xi, дискретно описывающих объем V(x):

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

Литература Гандин Л. С., Каган Р. Л. Статистические методы интерполяции метеорологических данных. — Л.: Гидрометеоиздат, 1976. — 359 с.

Armstrong M. Problems with universal kriging // Mathematical Geology. — 1984. — Vol. 16. — P. 101—108.

Chiles J. P., Delfiner P. Geostatistics. Modeling Spatial Uncertainty. — New York: A Wiley-Interscience Publication, 1999.

Dowd P. A. Lognormal kriging — the general case // Mathematical Geology. — 1982. — Vol. 14. — P. 475—499.

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

Oxford Univ. Press, 1997. — 483 p.

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

Kanevski M. F., Arutyunyan R. V., Bolshov L. A. et al. Spatial data analysis of Chernobyl fallout data. — 1. Preliminary results / Nuclear Safety Inst. — Moscow, 1993. — 91 p. — (Preprint NSI-23-93).

Kanevski M., Arutyunyan R., Bolshov L. et al. Mapping of Radioactively Contaminated Territories with Geostatistics and Artificial Neural Networks // Contaminated Forests / Eds. I. Linkov and W. R. Schell. — [S. l.]: Kluwer Academic Publ., 1999. — P. 249—256.

Krige D. G. A statistical approach to some basic mine valuation problems on the Witwatersrand // J. of the Chem., Metal. and Mining Soc. of South Africa. — 1951. — Vol. 52. — Р. 119—139.

Глава Многопеременное пространственное моделирование Данная глава посвящена проблемам использования дополнительной информации в рамках классической геостатистики. В Разделе 6.1 рассмотрены особенности использования информации, известной на всей области (кригинг с внешним дрейфом). В Разделах 6.2, 6.3 описаны взаимная пространственная корреляция нескольких переменных и их совместное оценивание (кокригинг). Уменьшение размерности пространства переменных с помощью метода принципиальных компонент и переход к факторному кригингу рассмотрены в Разделе 6.4.

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

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

6.1. Кригинг с внешним дрейфом Кригинг с внешним дрейфом (kriging with external drift) можно рассматривать как модификацию универсального кригинга (см. Раздел 5.4). В универсальном кригинге тренд моделируется линейной комбинацией базисГлава ных функций (5.20), в то время как в кригинге с внешним дрейфом тренд моделируется как линейная функция гладкой дополнительной переменной y(x), внешней по отношению к оцениваемой Z(x):

Неизвестные коэффициенты предполагаются постоянными в окрестности оцениваемой точки и вычисляются при решении кригинговой системы уравнений. В оценку Z*(x) внешний дрейф входит опосредованно через влияние на весовые коэффициенты j(x), j = 1,..., n(u):

Система уравнений кригинга с внешним дрейфом при этом имеет вид где CR(xj – x) — ковариационная функция невязки R(x) = Z(x) – m(x). Система уравнений кригинга с внешним дрейфом является частным случаем системы уравнений универсального кригинга (5.21), если K = 1 и f1(x) в любой точке совпадает со значением вторичной переменной y(x).

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

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

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

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

В. В. Демьянов, Е. А. Савельева Геостатистика: теория и практика • Должна быть возможность оценки и моделирования ковариации CR(xj – x) или вариограммы R(xj – x) остатков по значениям реальных измерений, так как сами остатки становятся известны только после решения кригинговой системы. Это требование непосредственно связано с предыдущим о плавности изменений вторичной переменной — y(xj) y(xj + h).

Примером использования кригинга с внешним дрейфом может служить моделирование поля температуры при наличии дополнительной информации о высоте над уровнем моря на подробной сетке (практически в любой точке). Такая информация доступна в виде цифровой карты уровней (digital elevation map). На рис. 6.1 приведены карта пространственного распределения данных о температуре и цифровая карта уровней для этого места. Рисунок 6.2 показывает практически линейную зависимость между температурой и высотой. Такая ярко выраженная линейная зависимость позволяет использовать данные о высоте как внешний дрейф для температуры. Применение обычного кригинга на специально выбранном из исходных данных валидационном наборе дало среднеквадратичную ошибку 3,13. Использование дополнительной информации позволило уменьшить среднеквадратичную ошибку на валидационном наборе до 1,42.

Рис. 6.1. Точки измерений температуры на поверхности (а), Рис. 6.2. Зависимость между температурой воздуха у поверхности 6.2. Меры корреляции и пространственной корреляции нескольких переменных Теперь рассмотрим общий случай измерений нескольких параметров, интерпретируя их как многопеременную случайную функцию. Итак, пусть Z(xi) — многопеременная функция, заданная на области S (i = 1,..., n; = 1,..., K), где i — индекс, означающий номер измерения (пространственной точки);

— индекс, означающий номер переменной. Все измерения всех переменных можно представить в виде матрицы Z размерностью Kn.

Корреляция переменных описывается ковариационной матрицей:

где m — вектор средних значений отдельных переменных.

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

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

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

Кросс-ковариация (cross-covariance). Для N(h) экспериментальных точек, разделенных вектором h, в которых есть измерения обеих переменных Z(x) и Z(x), кросс-ковариация определяется как [Dowd, 1989] где mZ h — среднее значение переменной Z по началу вектора h, а mZ +h — среднее значение переменной Z по концам вектора h.

Кросс-ковариационная функция не является априори ни симметричной, ни инвариантной относительно перестановки переменных:

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

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

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

Кросс-вариограмма (cross-variogram). Кросс-вариограмма для N(h) экспериментальных точек, разделенных вектором h, в которых есть измерения обеих переменных Z(x) и Z(x), определяется по формуле Как и обычная вариограмма, кросс-вариограмма обладает симметрией, кроме того, она инвариантна относительно перестановки переменных:

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

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

• полная гетеротопия (complete heterotopy) — измерения переменных находятся в различных точках и не имеют ни одной общей точки;

• частичная гетеротопия (partial heterotopy) — переменные имеют и общие, и различные точки измерений;

• изотопия (isotopy) — в каждой точке сети мониторинга есть данные измерений по всем переменным.

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

Другим путем при полной гетеротопии может быть использование псевдокросс-вариограммы, которая учитывает значения дополнительной переменной в точках, в которых отсутствуют данные по основной переменной [Papritz et al., 1993]. Для учета всех имеющихся данных в совпадающих и не совпадающих точках измерений переменных даже в случае полной гетеротопии псевдокросс-вариограмму можно вычислить следующим образом:

где xi, yi — векторы координат данных, разделенных вектором h (см. Раздел 4.3).

Псевдокросс-вариограмма может быть выражена через обыкновенную кросс-вариограмму:

Значения псевдокросс-вариограммы g могут быть использованы в уравнениях кокригинга (см. Раздел 5.3), как и обыкновенная кросс-вариограмма.

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

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

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

Линейная модель корегионализации предполагает, что каждая компонента вектора случайной функции может быть представлена как линейная комбинация некоррелированных компонент. Они обычно представляются в виде моделей ковариации или вариограмм одного типа, но с разными радиусами корреляции. Преимущество линейной модели корегионализации состоит в том, что условие положительной определенности сводится к проверке положительной определенности постоянной матрицы. Она наиболее полезна при малом количестве измерений. Ее недостаток — ограниченный выбор моделей кросс-вариограмм и кросс-ковариаций. Другими словами, если каждая компонента представлена как линейная комбинация некоррелированных компонент, то это соответствует диагонализации матрицы структурной функции [Myers, 1995].

Рассмотрим вектор значений случайной функции Z(x) = [Z1(x),..., Zm(x)].

Пусть Y1(x),..., Yp(x) — некоррелированные случайные функции, где p может быть больше или меньше m. Стационарность компонент Y следует из стационарности компонент Z и наоборот. Предположим, что Выражение (6.5) можно представить в матричной форме:

Пространственные функции Z соотносятся с пространственными функциями Y следующим образом:

где CZ(h), C(h) — ковариационные матрицы с компонентами:

а Z(h), Y(h) — матрицы вариограмм с компонентами:

uv,Y(h) = Cov{Yu(x + h) – Yu(x), Yv(x + h) – Yv(x)} = 0, при u v. (6.10) Заметим, что в уравнениях для вариограмм (6.9) и (6.10) Z(h) и Y(h) являются стандартными кросс-вариограммами, а не приведенными выше псевдокросс-вариограммами — см. (6.4) [Myers, 1991].

Из (6.7) и (6.8) для ковариации получим или Линейная модель корегионализации обычно записывается в виде (6.12).

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

Для случая двух переменных U и V модели авто- и кросс-вариограмм строятся следующим образом:

В. В. Демьянов, Е. А. Савельева Геостатистика: теория и практика где U(h), (h), UV(h) — авто- и кросс-вариограммные модели для U и V соответственно. Базисные модели задаются 0(h), 1(h),..., m(h); u, v, w — коэффициенты, возможно отрицательные.

Уравнения (6.13) можно записать в матричной форме:

комбинация первых базисных моделей 0(h):

комбинация вторых базисных моделей 1(h):

комбинация m-х базисных моделей m(h):

Для удовлетворения условия положительной определенности линейной модели (6.13) достаточно положительности коэффициентов u, v, w в уравнениях (6.14)—(6.16). Это достигается наложением на коэффициенты следующих условий:

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

Рис. 6.3. Экспериментальная и модельные вариограммы 90Sr (а), 137Cs (б) и кроссвариограмма 90Sr-137Cs (в) с применением линейной модели корегионализации В. В. Демьянов, Е. А. Савельева Геостатистика: теория и практика Рисунок 6.3 и табл. 6.1 иллюстрируют разрешение проблемы корегионализации при совместном моделировании пространственной структуры 137Cs и 90Sr.

Для всех моделей выполняется условие положительной определенности:

для наггета:

для параметров сферической модели:

Таблица 6.1. Сферические модели авто- и кросс-вариограмм для 90Sr и 137Cs 6.4. Кокригинг Кокригинг (cokriging) — естественное обобщение кригинга на случай многопеременных данных, когда между переменными имеется пространственная корреляция. Основная переменная оценивается на основе ее собственных измерений и данных по другим (дополнительным) переменным. Знание всех переменных во всех точках не требуется. Для обычного кокригинга обязательно по крайней мере одно измерение основной переменной, для простого же достаточно знания ее среднего значения, остальная информация вносится за счет дополнительных переменных. С другой стороны, случай полной изотопии данных даже при взаимной коррелированности переменных эквивалентен кригингу и не дает дополнительного улучшения оценки.

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

Оценка обычного кокригинга функции Z ( x0 ) — линейная комбинация значений различных переменных из окрестности точки x0. Количество участников оценивания среди различных переменных n может быть различно:

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

Рассмотрим условие несмещенности для так построенной оценки:

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

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

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

Она может быть выражена в терминах кросс-ковариации и ковариации или в терминах кросс-вариограммы и вариограммы. Здесь приведен вариант выражения вариограммы, который более характерен при использовании обычного кокригинга (когда неизвестно среднее) [Wackernagel, 1995]:

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

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

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

где m — стационарные средние значения Z(x). Для реализации условия несмещенности (2.12) в этом случае используется другое условие:

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

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

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

Пример использования кокригинга для Чернобыльских данных. Полезность кокригинга удобно проиллюстрировать на примере анализа пространственных данных по загрязнению окружающей среды в результате ЧерноГлава быльской аварии изотопами 137Cs и 90Sr, которые сильно коррелированны между собой (рис. 6.4). Если рассматривать их корреляцию как линейную, то коэффициент корреляции равен 0,74. Количество измерений 90Sr было меньше, чем 137Cs, из-за их дороговизны (рис. 6.5). Методика кокригинга позволила использовать измерения 137Cs для уточнения оценки загрязнения Sr. В частности, зоны, где отсутствуют измерения по 90Sr, но присутствуют измерения по 137Cs (северо-восточная часть данной территории), перестают быть зонами экстраполяции. Именно для них использование кокригинга и представляет особый интерес, так как там, где есть измерения по 90Sr, достаточно оценки обычного кригинга.

Рис. 6.4. Корреляция между значениями выпадений 137Cs и 90Sr При использовании обычного кригинга оценивание в областях экстраполяции связано с выбором области поиска используемых при оценке данных. Размер анизотропной корреляционной структуры для 90Sr (до 40 км, см. вариограммную поверхность на рис. 6.6б) не позволяет получить корректную оценку в этих областях экстраполяции. Кокригинг же использует для оценивания дополнительные измерения 137Cs, присутствующие в этих областях. Кроме того, более протяженная корреляционная структура 137Cs (до 70 км, см. вариограммную поверхность на рис. 6.6а) и соответствующая ей область поиска дают возможность захватить более обширную территорию. Таким образом, кокригинг позволяет избежать появления не оцененных (или оцененных некорректно) областей по сравнению с обычным кригингом.

На рис. 6.7 представлена пространственная кросс-корреляционная структура 137Cs и 90Sr. Пространственная корреляция и пространственная кросскорреляция представлены в виде вариограммных (кросс-вариограммных) В. В. Демьянов, Е. А. Савельева Геостатистика: теория и практика поверхностей. На рисунках 6.6, 6.7 приведены и экспериментальные, и соответствующие им модельные структуры.

Рис. 6.5. Схема точек снятия проб по загрязнению поверхности Рис. 6.6. Экспериментальные (верхняя и нижняя) и модельные (средняя) вариограммные поверхности для данных по загрязнению поверхности западной части Брянской области:

Рис. 6.7. Экспериментальные (верхняя и нижняя) и модельная (средняя) кросс-вариограммные поверхности для данных по загрязнению 137Cs и 90Sr поверхности западной части Брянской области Кокригинг выполнялся на прямоугольной сетке 4570 ячеек с размером ячейки 22 км. Полученная карта оценок 90Sr представлена на рис. 6.8. Для сравнения оценка загрязнения 90Sr почвы в этой же области с использованием обычного кригинга представлена на рис. 6.9. При сравнении карт оценок видно, что обычный кригинг сделал размазывающее усреднение к границам и оставил не оцененными значительные области по краям. Это связано с отсутствием в этих областях измерений 90Sr. Кокригинг же использует для оценивания дополнительные измерения 137Cs в областях, где измерения 90Sr отсутствуют.

Рис. 6.8. Карта результатов оценки загрязнения 90Sr западной части Брянской области с помощью обычного кокригинга В. В. Демьянов, Е. А. Савельева Геостатистика: теория и практика Рис. 6.9. Карта результатов оценки загрязнения 90Sr западной части Брянской области с помощью обычного кригинга Оценка обычного кригинга более размазана, менее вариабельна. Это означает занижение оценки в точках с высоким загрязнением и ее завышение в точках с низким загрязнением. В отличие от обычного кригинга кокригинг дал значительно менее сглаженную оценку. На карте четко видны характерные пятна большого загрязнения вблизи точек выбросов высоких значений.

Сравнение итоговых статистик оценки кокригинга и исходных измерений (с учетом и без учета декластеризации, см. Раздел 2.5) также подтверждает ее превосходство над оценкой обычного кригинга (табл. 6.2). Отметим, что кокригинг позволяет улучшить оценку высоких значений.

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

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

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

Таблица 6.2. Статистика распределений измерений Количество данных Среднее значение Вариация Стандартное отклонение Коэффициент вариации, % Минимум Нижний квартиль (25%) Медиана Верхний квартиль (75%) Максимум Пример использования кокригинга для исследования загрязнения Женевского озера. Можно сравнить результаты применения кригинга и кокригинга с различным количеством переменных для случая девятипеременной функции (загрязнение донных отложений Женевского озера металлами).

Одна из переменных является основной (Pb), остальные (Zn, Cu, Mn, Cd, Ni, Be, B, Cr) — дополнительной информацией. На рис. 6.10 представлены кросс-вариограммы основной переменной со всеми дополнительными переменными. Видно, что пространственная корреляция для всех, кроме одной (правая верхняя с B), достаточно хорошо моделируется.

Рис. 6.10. Кросс-вариограммы основной переменной (Pb) с дополнительными переменными. Линии, соединяющие значки, — экспериментальные кроссвариограммы, кривые — модели кросс-вариограмм В. В. Демьянов, Е. А. Савельева Геостатистика: теория и практика Сравнение результатов кокригинга проводилось по среднеквадратичной ошибке кросс-валидации. Результаты собраны в табл. 6.3. Все кокригинги дают меньшую ошибку кросс-валидации, чем кригинг (86,37). Использование третьей переменной улучшает оценку по сравнению с двумя. Но введение четвертой переменной (при использовании различных наборов) уже ведет к ухудшению. Использование всех переменных дает худший для кокригинга результат.

Таблица 6.3. Среднеквадратичные ошибки кросс-валидации (RMSE) при оценке кокригинга основной переменной (Pb) с использованием различных наборов дополнительных переменных (Zn, Cu, B, Cd, Cr, Mn, Be, Ni) Количество дополнительных Другой причиной того, что кокригинг не очень популярен, является эффект экранирования более коррелированными данными (измерениями основной переменной) менее коррелированных (измерений дополнительной переменной). Влияние эффекта экранирования на значения весов и их отрицательный эффект на оценку уже обсуждались при рассмотрении обычного кригинга (см. Подраздел 5.6.1). В случае кокригинга эффект экранирования проявляется чаще за счет дополнительной информации в дополнительных точках измерения.

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

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

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

В Разделе 6.1 уже упоминался пример моделирования поля температур с использованием дополнительной информации о высоте над уровнем моря.

Применение полного кокригинга в этом примере привело к среднеквадратичной ошибке на валидационном наборе в размере 3,97. Это даже хуже, чем результат обычного кригинга, где среднеквадратичная ошибка на валидационном наборе составляла 3,13. Использование колокационного кокригинга позволило уменьшить среднеквадратичную ошибку на валидационном наборе до 3,05.

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

Одним из способов уменьшения вычислительной сложности многопеременного кокригинга является использование кригинга принципиальных компоВ. В. Демьянов, Е. А. Савельева Геостатистика: теория и практика нент. Этот подход предполагает переход от K коррелированных переменных к d (d K) некоррелированным факторам. Эти факторы представляют собой линейную комбинацию исходных переменных и могут быть использованы для описания данных при меньшей размерности. Для перехода к факторам используется известный в многопеременной статистике анализ принципиальных компонент [Вентцель, 1964], где эти принципиальные компоненты и являются искомыми факторами.

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

Рис. 6.11. Схема построения новой системы координат Рассмотрим анализ принципиальных компонент в математическом плане. Пусть имеется матрица данных, из которых вычтены значения средних по каждой переменной. Обозначим ее через Z, она имеет размерность KN — N точек измерений по K переменных в каждой. Для нее можно построить ковариационную матрицу VZ по формуле (6.1) (среднее равно нулю). Требуется найти такую ортогональную матрицу A размерности NN, что в результате ее умножения на исходную матрицу Z получим новую матрицу Y, обладающую диагональной ковариационной матрицей VY, т. е.

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

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

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

На рис. 6.12 показаны значения принципиальных компонент для примера, уже использовавшегося в этой главе (содержание металлов в донных отложениях Женевского озера). Пунктирная линия проводит границу значимости (90%). Три компоненты вносят основной вклад. Этот результат согласуется с тем, что был получен при использовании кокригинга, — три переменные (если они были правильно выбраны) давали лучший результат.

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

Рис. 6.13. Круги корреляции девяти переменных для плоскостей, определяемых факторами 1 и 2 (слева) и факторами 1 и 3 (справа) Для каждого из значимых факторов, полученных в результате анализа принципиальных компонент, может быть проведен кригинг. Так как факторы по определению ортогональны (и в соответствии с нашим желанием не коррелированны), их совместную оценку проводить бессмысленно. По оцененным факторам можно вернуться к исходным данным, используя обратное преобразование Z = A–1Y. Таким образом, вместо моделирования K моделей пространственной корреляции и кросс-корреляции мы обходимся моделированием d K пространственных корреляций значимых факторов, но в дополнение решаем задачу поиска собственных векторов и собственных значений ковариационной матрицы.

Литература Вентцель Е. С. Теория вероятностей. — М., 1964.

Каневский М. Ф., Арутюнян Р. В., Большов Л. А. и др. Геостатистический подход к анализу чернобыльских выпадений // Изв. Рос. акад. наук.

Энергетика. — 1995. — № 3. — С. 34—46.

Cressie N. Statistics for Spatial Data. — New York: John Wiley & Sons, 1991. — Р. 141.

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

Dowd P. A. Generalised cross-covariances // Geostatistics. — 1989. — Vol. 1. — P. 579—590.

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

Myers D. E. Pseudo-Cross Variograms, Positive-Definiteness and cokriging // Mathematical Geology. — 1991. — Vol. 23, N 6. — P. 805—816.

Myers D. E. The Linear coregionalization and simultaneous diagonalization of the variogram matrix function // Sciences de la Terre. — 1995. — Vol. 32. — Р. 125—139.

Pan G., Gaard D., Moss K., Heiner T. A Comparison Between Cokriging and Ordinary Kriging: Case Study with a Pollymetallic Deposit // Mathematical Geology. — 1993. — Vol. 25, N 3. — P. 377—398.

Papritz A., Kunsch H. R., Webster R. On the Pseudo Cross Variogram // Mathematical Geology. — 1993. — Vol. 25, N 8. — Р. 1015—1026.

Wackernagel H. Multivariate Geostatistics. — Berlin: Springer-Verl., 1995.

Глава Вероятностное моделирование локальной неопределенности Эта глава включает описание индикаторного преобразования (Раздел 7.1), применения индикаторного подхода для анализа непрерывной и категориальной переменных (Раздел 7.2). Рассмотрены примеры использования индикаторного подхода для различных типов данных (Раздел 7.3).

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

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

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

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

Индикаторный подход был изложен в работах [Journel, 1983, 1985]. Было предложено использовать при анализе и моделировании переход от исходных значений переменных к специальным индикаторам. Индикаторные переменные — бинарные, т. е. принимают значения либо 0, либо 1. Для категориальной переменной такое преобразование дает индикаторную Вероятностное моделирование локальной неопределенности переменную для каждой из категорий, характеризуя присутствие или отсутствие данной категории. В случае непрерывной переменной переход представляет собой нелинейное преобразование, моделирующее кумулятивную функцию распределения. Переход к индикаторным переменным позволяет получить оценку условной локальной функции распределения (но не ее явный вид) во всем пространстве. Преобразованные данные более устойчивы к выбросам (outliers). Индикаторное преобразование может быть применено и при анализе категориальных данных, а также при совместном анализе данных различных типов, что, в свою очередь, добавляет такому подходу дополнительную привлекательность.

7.1. Индикаторное преобразование Для индикаторного преобразования непрерывной случайной функции Z(x) сначала выбирается набор пороговых значений zk, k = 1,..., K. Затем производится переход к индикаторам, для каждого порогового значения zk определяется индикатор В результате для каждой точки x пространства получается вектор индикаторов (размерности K):

где K — число пороговых отсечений.

Пример проведения индикаторного преобразования, определяющего тип фации (породы) для одного порогового значения, представлен на рис. 7. на примере профиля каротажа Z в скважине. Разбиение на категории осуществляется в соответствии с пороговым значением zk — его превышением и непревышением. Значения функции больше порогового значения трансформируются в значения 0, значения функции меньше порогового значения преобразуются в 1. Использование набора пороговых значений в случае непрерывной переменной дает возможность построить подробную локальную функцию распределения с требуемым разрешением.

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

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

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

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

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

Вероятностное моделирование локальной неопределенности Теперь рассмотрим функцию Q(x), определенную на области S и принимающую конечный набор значений (c1,..., cC). Такую функцию называют пространственной (региональной) категориальной. На практике это могут быть типы почв, типы геологического формирования, предупредительный сигнал прибора и т. п. Для пространственной категориальной переменной Q(x) индикаторное преобразование делается для каждого возможного значения (класса) ci, i = 1, …, C:

Для категориальной переменной вектор индикаторов состоит из 0 и одной 1 в соответствии со значением Q(x). Глобальное усреднение индикаторов даст относительное распределение исходного набора данных по классам (значениям категориальной функции):

Индикаторы и непрерывной, и категориальной переменных по сути относятся к одному типу — это категориальные индикаторы, принимающие два значения — 0 и 1. Поэтому их можно использовать вместе.

Отдельный индикатор — I(x, zk) или I(x, ci) — можно рассматривать как пространственную функцию от x. Соответственно можно ввести и пространственные корреляционные функции для этих переменных.

Нецентральная индикаторная ковариация:

В. В. Демьянов, Е. А. Савельева Геостатистика: теория и практика Центральная индикаторная ковариация:

• для индикатора непрерывной переменной:

• для индикатора категориальной переменной:

Индикаторная полувариограмма:

• для индикатора непрерывной переменной:

• для индикатора категориальной переменной:

На практике при использовании индикаторного подхода для анализа пространственных данных используется индикаторная полувариограмма. Анализ и моделирование индикаторной полувариограммы проводится по методике, описанной в Главе 4.

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

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

В Главе 5 рассматривался пример вариограммы для траловой съемки краба Берди в Беринговом море в 2003 г. Там разброс значений составлял от 0 до 170000, (т. е. пять порядков), что не позволяло оценивать пространственную корреляцию с помощью вариограммы (см. рис. 5.16). ИндикаГлава Вероятностное моделирование локальной неопределенности торные вариограммы для четырех порогов, соответствующих квантилям 0, (321 краб), 0,625 (870,5 краба), 0,75 (2186 крабов) и 0,875 (7421,13 краба), представлены на рис. 7.2. Каждая из них имеет структуру, которую вполне можно моделировать.

Рис. 7.2. Индикаторные вариограммы для данных траловой съемки пространственного распределения краба Берди 7.2. Индикаторный кригинг Индикаторным кригингом (indicator kriging) называется обычный кригинг, выполненный для индикаторов, т. е. это линейный оцениватель, построенный по аналогии с обычным кригингом, но не для значений анализируемой переменной, а для индикатора:

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

Веса ki, полученные при решении системы уравнений обычного кригинга для индикаторов, В. В. Демьянов, Е. А. Савельева Геостатистика: теория и практика удовлетворяют требованию минимизации вариации ошибки. Здесь n — число точек измерений с заданными значениями оцениваемой функции Z(xi); x0 — точка, в которой производится оценивание; hij = xj – xi — вектор разности координат между соответствующими точками, который является аргументом функций ковариаций и вариограмм (по условию стационарности для Z(x)).

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

Рис. 7.3. Пример вариограмм, связанных множителями с базовой вариограммой Здесь выделен базовый порог zM, через модель пространственной корреляции которого описываются пространственные корреляции всех остальных отсечений. Как уже упоминалось (см. Главу 5), веса кригинга i не изменяются при мультипликативном преобразовании модели пространственной корреляции. Поэтому, будучи вычислены один раз, они могут быть использованы и для получения оценки индикатора, соответствующего другому порогу. Этот упрощенный вариант индикаторного кригинга принято называть медианным кригингом, так как обычно при выполнении условия (7.6) в качестве базового порога выбирается значение медианы исходного набора данных.

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

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

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

• сначала делается коррекция на попадание в отрезок [0, 1]:

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

• при проходе вниз проводится коррекция в сторону уменьшения, если наблюдается рост оцененного значения при уменьшении порогового значения:

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

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

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

где — положительная величина. При = 1 получается линейная интерполяция. При < 1 более быстрый, чем линейный, рост оценки функции распределения в начале интервала сменяется более медленным на его второй половине. Чем меньше значение, тем дальше оценка функции распределения от линейной. При > 1 все происходит наоборот. Примеры такой степенной оценки для одного отрезка интерполяции приведены на рис. 7.4.

Рис. 7.4. Примеры степенной интерполяции, применяемой Вероятностное моделирование локальной неопределенности При необходимости оценить значение функции распределения за пределами минимального или максимального порогов (случай экстраполяции) можно пользоваться тем же методом, но требуется задать локальные значения zmin и zmax. Для случая оценки меньше минимального порога (нижний хвост) используется формула а при z больше максимального порога (верхний хвост) — К этим случаям относится все, что было сказано выше относительно параметра степенной экстраполяции (см. рис. 7.4).

Для оценки верхнего хвоста можно использовать гиперболическую модель экстраполяции. Это позволяет строить функцию распределения до бесконечности, асимптотически приближая ее к единице. Гиперболическая модель имеет вид где параметр (> 1) контролирует скорость приближения к граничному значению функции распределения (т. е. к единице). Чем больше, тем короче хвост. Параметр определяется оценкой функции распределения, полученной для максимального порога:

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

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

Оценки Е-типа (E-type) — среднее значение в точке. Название происходит от английской транскрипции математического ожидания (expectation). Они вычисляются по формуле Эти оценки можно сравнивать с интерполяциями другими методами, например обычным кригингом. Средние оценки Е-типа можно представить на карте и гистограмме. Как всякое усреднение, оценка Е-типа дает сглаженное значение.

P-квантиль — вероятность превышения (или непревышения 1 – p) некоторого уровня значений zc, что соответствует оцененной функции распределения. Понять это можно следующим образом: квантиль p = 0,1 значит, что настоящее (но неизвестное) значение функции с вероятностью 90% превысит этот уровень. Одновременно можно рассчитать средние значения оценок выше и ниже порога отсечения. Картирование таких вероятностей может быть полезно при анализе последствий выбросов для поддержки принятия квалифицированных решений, например о проведении контрмер.

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

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

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

Коррекция для удовлетворения равенства суммы единице состоит в традиционной нормировке:

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

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

• Системы уравнений кокригинга также имеют большую размерность.

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

Хотя в теории индикаторный кокригинг должен быть лучше, на практике это не подтверждается [Goovaerts, 1997].

7.3. Примеры использования индикаторного подхода 7.3.1. Зонирование гидрогеологического слоя В этом примере рассматривается применение индикаторного кригинга к категориальным данным — зонирование гидрогеологического слоя. Зонирование — это подход, альтернативный анализу пространственной вариабельности гидрологических параметров, таких как пористость, проводимость и т. п. Он состоит в том, что область разбивается на зоны, гидрологические параметры внутри которых считаются постоянными. Обычно использование такого подхода обусловлено недостатком измерений для моделирования непосредственно гидрологических параметров.

В рассматриваемом примере проводится зонирование хорошо проводящего гидрогеологического слоя [Savelieva et al., 2003], который представлен пятью типами (зонами): тремя типами гравия, песка и ила. Исходные данные представляют собой 225 измерений, которые изображены на рис. 7.6 с помощью полигонов Вороного, которые здесь использованы исключительно для улучшения визуализации, так как при рисовании пространственного распределения точек измерений часть из них перекрывается.

Вероятностное моделирование локальной неопределенности Рис. 7.6. Исходные данные по гидрогеологическим классам Первый шаг состоит в индикаторном преобразовании исходных данных:

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

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

В. В. Демьянов, Е. А. Савельева Геостатистика: теория и практика Окончательное решение (классификация) принимается в пользу класса с максимальной вероятностью. На рис. 7.8 приведены результат классификации и вероятность класса-победителя, т. е. максимальная среди классов в этой точке.

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

На рис. 7.9 представлен результат классификации для задачи зонирования гидрогеологического слоя с использованием такого правила классификации и границей допустимой вероятности классификации 0,7. На рисунке видны белые зоны — это неклассифицированные зоны. Видно, что все зоны неопределенности расположены вдоль границы смены классов, следовательно, они являются аналогом «толстых изолиний» (подробнее о «толстых изолиниях» см. в Главе 5); где-то в пределах этих зон происходит смена одного класса на другой.

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

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

Рис. 7.10. Пространственное распределение индикатора для порогового значения 5000 крабов опилио Индикаторный кригинг позволяет оценить вероятность найти скопление крабов, т. е. вероятность, что число крабов в данном месте будет больше или равно 5000. Оценка индикаторного кригинга дает вероятность того, что значение функции не превышает порогового значения. Но если вычесть оценку индикаторного кригинга из единицы, то получится искомая величина. На рис. 7.11 изображена вероятность обнаружения скопления краба опилио. Белые плюсы показывают места, где результаты измерений указывали на скопления. Соответствие оценки и реальных измерений представляется вполне хорошим.

В. В. Демьянов, Е. А. Савельева Геостатистика: теория и практика Рис. 7.11. Карта вероятности обнаружить скопление крабов опилио Литература Goovaerts P. Geostatistics for Natural Resources Evaluation. — [S. l.]:

Oxford Univ. Press, 1997. — 483 p.

Journel A. G. Nonparametristic Estimation of Spatial Distribution // Mathematical Geology. — 1983. — Vol. 15, N 3.

Journel A. G. The Deterministic Side of Geostatistics // Mathematical Geology. — 1985. — Vol. 17, N 1.

Savelieva E., Bolshov L., Pozdnukhov A. et al. Zonation of Hanford Formation at the Hanford Site / Nuclear Safety Inst. RAS. — Moscow, 2003. — 55 p. — (Preprint IBRAE-2003-07).

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

Раздел 8.1 посвящен основам стохастического моделирования и его принципиальному отличию от регрессионного оценивания. Здесь же представлены основные подходы к стохастическому моделированию, его виды и кратко перечислены существующие алгоритмы. В разделе 8.2 описан ключевой для большинства стохастических моделей геостатистики принцип последовательного моделирования. Следующие разделы посвящены конкретным алгоритмам, которые базируются на последовательном принципе моделирования. В разделе 8.3 подробно описано последовательное гауссово моделирование, в разделе 8.4 — обрезанное гауссово моделирование, в разделе 8.5 — последовательное индикаторное моделирование, в разделе 8.6 — последовательное прямое моделирование. Принципиально другой алгоритм — моделирование отжига — представлен в разделе 8.7. Раздел 8.8 посвящен объектному подходу к стохастическому моделированию.

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

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

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

Для оценки вариабельности пространственной функции используют методы стохастического моделирования, которое в отличие от кригинга позволяет получить множество реализаций значений функции в точке оценивания для заданного набора данных и выбранных параметров модели [Journel, Huijbregts, 1978].

Проиллюстрируем сказанное на простом одномерном примере (рис. 8.1).

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

Принципиальное различие между реализациями состоит в разнообразии возможных значений вне точек измерений. Вследствие стохастической природы модели функция может принимать значения выше и ниже измеренного и воспроизводить разнообразие стохастической динамики изменения значений — вариабельности. Отметим, что вариабельность меняется в зависимости от близости данных измерений. Так, в интервале абсцисс {8, 20} и {59, 76} вариабельность значительно ниже, чем там, где отсутствуют измерения.

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

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

• плотность распределения;

• статистические характеристики исходных данных;

• пространственную корреляционную структуру.

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

Стохастическое моделирование может быть условным — зависимым от данных (см. рис. 8.1) — или безусловным, когда нет данных измерений.

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

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

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

Самый известный метод безусловных симуляций — метод «вращающихся лент»

(turning bands) [Mantoglou, Wilson, 1981]. Он позволяет построить безусловные реализации гауссова поля Y(x) с известной ковариационной функцией CY(h).

Основная идея метода состоит в построении одномерных реализаций вдоль 15 линий, различным образом ориентированных и разбивающих трехмерную область на примерно равные части. Каждый узел 3D-пространства, где производится моделирование, проектируется в некоторую точку на каждой линии.

Значение в узле задается суммой значений этих проекций.

В. В. Демьянов, Е. А. Савельева Геостатистика: теория и практика Моделирование одномерной реализации (вдоль линии) обычно выполняется с использованием спектрального подхода на основе быстрого преобразования Фурье [Christakos, 1992]. Функция ковариации CY(h), описывающая пространственную корреляцию моделируемого поля, при этом подходе заменяется соответствующей ей функцией спектральной плотности SY():

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

Существуют различные подходы к стохастическому пространственному моделированию [Chiles, Delfiner, 1999]. Один из них основан на последовательном принципе моделирования. Другой подход использует методы глобальной минимизации целевой функции.

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

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

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

Пространственные реализации, полученные в результате применения любого метода стохастического моделирования, могут использоваться как входные данные модели оценки функции распределения пространственной переменной. Это дает возможность оценить локальную вариабельность пространственного распределения. Возможные реализации геологической среды, полученные в результате стохастического моделирования, могут использоваться для оценки вероятности попадания загрязнения через грунтовые воды. Реализации стохастического моделирования для ошибок при прогнозировании потребления электроэнергии могут использоваться при оценке неопределенности (доверительного интервала) прогноза и т. п. По результату стохастического моделирования можно оценивать вероятность превышения или непревышения определенных значений как в отдельной точке (задача моделирования локальной функции распределения), так и в нескольких точках одновременно (рассматривается совместная функция распределения). На практике чаще моделируются и используются двухточечные совместные функции распределения. Стохастическое моделирование позволяет также получить оценки вероятности превышения заданного уровня (p-value) и оценки с заданным уровнем вероятности превышения. Такие оценки, а также анализ неопределенности пространственной оценки, крайне важны для поддержки принятия квалифицированных решений. При усреднении стохастических реализаций можно получить среднюю оценку (E-type), сравнимую с аналогичной оценкой для индикаторного кригинга (см. Раздел 7.2). Разница между стохастическими реализациями позволяет оценить вариацию и разброс значений локальных функций распределения.

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

В. В. Демьянов, Е. А. Савельева Геостатистика: теория и практика Пиксельные модели • моделируют значение в каждой отдельно взятой ячейке;

• позволяют легко интегрировать дополнительную локальную информацию;

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

Объектные модели используют реалистичные структуры (объекты различной геометрической формы), но имеют следующие недостатки:

• предположение о формах объектов требует априорных знаний о структуре системы;

• возникают неопределенности, связанные с выбором форм, размеров и местоположения объектов;

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

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

1. Последовательное гауссово моделирование (sequential Gaussian simulation) — моделирует непрерывные переменные (например, значение загрязнения, пористости породы, биомассу рыбы).

2. Обрезанное гауссово моделирование (truncated Gaussian simulation) — моделирует категориальные переменные (например, типы почв, геологических пород).

3. Последовательное индикаторное моделирование (sequential indicator simulations) — моделирует как категориальные переменные, так и непрерывные (типы пород, уровень загрязнения).

4. Прямое моделирование (direct simulations) — моделирует непрерывные переменные.

5. Моделирование отжига (simulated annealing) — моделирует категориальные и непрерывные переменные.

6. Объектное моделирование (object modelling) — моделирует категориальные переменные, использует заданный набор геометрических объектов.