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

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

Детерминистические методы интерполяции предполагают наличие заданной аналитической зависимости между значениями функции в пространстве. Эти методы популярны из-за простоты вычисления оценки по заданной параметрической формуле. Наиболее широко применяемые «формульные» зависимости:

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

При использовании детерминистических методов предполагается, что анализируемые данные описываются некоторой детерминистической функцией Z(x, ), определенной на исследуемой области S, где x S — координаты точки; — набор внутренних параметров модели. Задача состоит в том, чтобы, базируясь на известных данных (Zi = Z(xi) — значения, измеренные в точках xi S) и на другой контекстной информации об исследуемом явлении, подобрать набор параметров и построить функцию Z(x, ) для всей исследуемой области S. После этого значение в любой точке просто вычисляется по формуле.

Детерминистические интерполяторы могут быть глобальными (все точки с известными значениями используются при интерполяции) или локальными Детерминистические методы пространственной интерполяции (только часть значений в точках, ближайших к оцениваемой, используются для интерполяции). Глобальные интерполяторы делают искомую функцию более сглаженной. При использовании локальных методов окрестность, используемая для оценки, может задаваться различными способами. Может быть фиксировано число ближайших к оцениваемой точке соседей, использующихся при интерполяции: N(x) = N = const. Тогда размер зоны, влияющей на значение в точке x, зависит от локальной плотности точек измерения. Возможно, наоборот, задание размера области (например, область поиска D), точки из которой используются при оценивании значения в x.

В этом случае N(x) будет меняться и в области с редкими измерениями при малом значении D могут возникнуть неоцененные зоны (в D-окрестности точки x недостаточно измерений).

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

3.1. Линейные интерполяторы Линейные интерполяторы представляют искомую функцию в виде линейной комбинации известных значений:

где Z*(x0) — оцениваемое значение в точке x0; Z(xi) — известные значения в точках измерений xi; N(x0) — количество исходных точек, принимающих участие в оценке для координаты x0; wi(x0) — весовые коэффициенты.

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

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

В. В. Демьянов, Е. А. Савельева Геостатистика: теория и практика Широко используется линейный интерполятор с весовыми коэффициентами, обратно пропорциональными расстоянию до оцениваемой точки в степени. Весовые коэффициенты, определяются по формуле где — степень; hi = d i 2 + 2 ; di — расстояние между точками xi и x0; — сглаживающий параметр. Нижняя часть дроби в весовом коэффициенте вводится для выполнения условия несмещенности оценки E Z * ( x0 ) Z ( x0 ) = 0, которое соответствует условию на весовые коэффициенты wi = 1.

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

В качестве значения степени чаще всего используется значение 2. Такой вариант метода известен как метод обратных квадратов [Grimm, Lynch, 1991].

Рис. 3.1. Результат интерполяции методом обратных квадратов:

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

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

Таким образом, выбор подходящего радиуса поиска тесно связан с понятиями пространственной непрерывности и пространственной корреляции, которые подробно рассмотрены в Главе 4.

Для подбора оптимального значения радиуса поиска можно использовать один из алгоритмов проверки качества оценки — кросс-валидацию или метод складного ножа, которые были описаны в главе 2. При использовании кросс-валидации оптимальный радиус поиска определяется путем минимизации кросс-валидационной ошибки. На этом принципе основаны некоторые алгоритмы пространственного автокартирования [Kanevski et al., 1999].

На рис. 3.2 изображены интерполяционные оценки методом обратных расстояний в степени на основе трех точек: (2, 2), (5, 8) и (8, 5).

Упражнение 3.1. Какие степени — 1, 2 или 3 — использованы для получения интерполяционных оценок А, B и C?

В. В. Демьянов, Е. А. Савельева Геостатистика: теория и практика Рис. 3.2. Оценки методом обратных расстояний в степени для степеней 1, 2 и Если область поиска определяется радиусом, то он может быть в явном виде введен в весовые коэффициенты, что реализовано, например, в весовых коэффициентах Крессмана [Grimm, Lynch, 1991]. Этот метод по сути аналогичен методу обратных квадратов:

где Di — радиус влияния i-й точки. Радиус влияния может быть задан постоянным по всей исследуемой области, а может меняться в зависимости от локальной плотности точек измерения. При использовании весовых коэффициентов Крессмана их рекомендуется нормализовать так, чтобы выполнялось условие Весовые коэффициенты можно определять и совсем иначе, например основываясь на системе полигонов Вороного (подробнее о полигонах Вороного см. в Главе 2). Такой подход называется методом естественного соседа.

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

Детерминистические методы пространственной интерполяции где pi0 — площадь зоны, арендованной точкой x0 у точки xi.

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

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

3.2. Полиномиальные методы Полиномиальные интерполяторы [Goodin et al., 1979] представляют значение в точке в виде полинома от координат. В двумерном случае — для точки x с координатами (x, y) Z*(x, y) = Pn(x, y), где Pn — полином n-й степени.

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

• плоскость: P ( x, y ) = a + bx + cy;

• билинейно-седловой: P.5 ( x, y ) = a + bx + cy + dxy;

• квадратичный: P2 ( x, y ) = a + bx + cy + dxy + ex 2 + fy 2 ;

• кубический: P3 ( x, y ) = a + bx + cy + dxy + ex 2 + fy 2 + gx 2 y + hxy 2 + В. В. Демьянов, Е. А. Савельева Геостатистика: теория и практика Теоретически можно использовать и полиномы более высокого порядка.

Они определяются максимальной степенью при x, максимальной степенью при y и совместной максимальной степенью. Все промежуточные степени в полиноме будут присутствовать.

Задача полиномиальной интерполяции сводится к тому, чтобы определить неизвестные коэффициенты ai так, чтобы полиномы максимально хорошо соответствовали данным в заданных точках. Для этого находят минимум по всем коэффициентам (a, b, c, d и т. д.) функции 2, задающей интегральную ошибку интерполяции и определяемой следующим образом:

Минимизация состоит в решении системы линейных уравнений с числом неизвестных, равным числу уравнений. Число уравнений (неизвестных) зависит от выбранного полинома.

Любой глобальный полиномиальный метод, вообще говоря, не является интерполятором в строгом смысле, скорее он относится к аппроксиматорам.

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

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

3.3. Метод базисных функций Оценка методом базисных функций строится как линейная комбинация из базисных функций:

где h0i — расстояние от точки x0 до точки xi; B(h0i) — базисная функция, определяемая от расстояния; ci — весовые коэффициенты. Коэффициент ci определяет алгебраический знак вхождения соответствующего члена и степень его влияния. Классический вариант метода является точным, но возможно введение сглаживающего параметра.

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

• обратный мультиквадрик: B( h ) = ;

• мультилогарифмический: B (h ) = lg h 2 + 2 ;

• мультиквадратичный: B( h ) = h + ;

• тонкий сплайн: B( h ) = ( h + ) lg( h + ).

В. В. Демьянов, Е. А. Савельева Геостатистика: теория и практика Весовые коэффициенты ci определяются из условия точности оценки в известных точках, т. е. во всех заданных точках (xi, yi) модель интерполяции должна давать оценку, равную заданным значениям V(xi, yi). Определение весов производится при = 0, при использовании сглаженного варианта используется при оценке. Таким образом, чтобы найти коэффициенты ci, необходимо решить систему из N линейных уравнений с N неизвестными:

Метод базисных функций обладает универсальностью и эффективностью.

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

Рис. 3.5. Результат интерполяции методом базисных функций Детерминистические методы пространственной интерполяции Литература Bartier P. M., Keller C. P. Multivariate interpolation to incorporate thetic surface data using inverse distance weighting // Computers and Geosciences. — 1996. — Vol. 22, N 7. — Р. 795—799.

Franke R. Scattered Data Interpolation: Test of Some Methods // Mathematics of Computation. — 1982. — Vol. 38, N 157. — Р. 181—200.

Goodin W. R., McRae G. J., Seinfeld J. H. A comparison of interpolation methods for sparse data: application to wind and concentration fields // J. of Applied Meteorology. — 1979. — Vol. 18. — P. 761—771.

Grimm J. W., Lynch J. A. Statistical analysis of error in estimating wet deposition using five surface estimation algorithms // Atmospheric Environment. — 1991. — Vol. 25A. — P. 317—127.

Kanevski M., Demyanov V., Chernov S. et al. Geostat Office for Environmental and Pollution Spatial Data Analysis // Mathematische Geologie. — 1999. — Vol. 3, N 4. — Р. 73—83.

Lattuada R., Raper J. Applications of 3D Delaunay triangulation algorithms in geoscientific modelling // http://www.iah.bbscr.ac.uk/phd/gisruk95.html.

Macidonio G., Pareschi M. T. An algorithm for the triangulation of arbitrarily distributed points: applications to the volume estimate and terrain fitting // Computers and Geosciences. — 1991. — Vol. 17, N 7. — Р. 859—874.

Okabe A., Boot B., Sugihara K. Spatial Tessellations: Concepts and Applications of Voronoi Diagrams. — New York: J. Wiley & Sons, 1992. — 532 p.

Saunderson H. C. Multiquadric surfaces in C // Computers and Geosciences. — 1994. — Vol. 20, N 7/8. — Р. 1103—1122.

SURFER (R) Version 8.0: Surface Mapping System / Golden Software, Inc. — [S. l.], 2002.

Tabois G. Q., Salas J. D. A Comparative Analysis, for Spatial Interpolation of Precipitation // Water Resources Bul. — 1985. — Vol. 21, N 3. — Р. 365 —380.

Weber D., Englund E. Evaluation and Comparison of Spatial Interpolators // Mathematical Geology. —1992. — Vol. 24, N 4. — Р. 381—391.

Глава Анализ и моделирование пространственной корреляции. Вариография В этой главе подробно рассмотрена основная тема геостатистики — вариография. Под вариографией понимают анализ и моделирование пространственной корреляционной структуры данных. В Разделе 4.1 мы еще раз вернемся к пространственной непрерывности, описанной в Разделе 2.6.

В Разделе 4.2 собраны различные меры пространственной корреляции, количественно характеризующие пространственную непрерывность. Раздел 4. посвящен построению экспериментальной вариограммы для набора пространственных данных. Моделирование построенной экспериментальной вариограммы с помощью аналитических функций описано в Разделе 4.4.

В Разделах 4.5, 4.6 изложены свойства вариограммы на больших расстояниях и вблизи нуля. Различные типы анизотропной пространственной корреляции описаны в Разделе 4.7. В Раздел 4.8 вынесено практическое упражнение на определение вариограмм для различных пространственных образов. Проблемы, связанные с использованием вариограммы, рассмотрены в Разделах 4.8 и 4.9. В Разделе 4.10 приведен пример анализа и моделирования пространственной корреляционной структуры для реальных данных.

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

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

Анализ пространственной корреляционной структуры данных можно разбить на два этапа:

• построение и интерпретация мер пространственной непрерывности на основе данных;

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

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

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

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

Приняв предположение о стационарности второго порядка или внутреннюю гипотезу, считаем, что функции корреляции между данными зависят только от взаимного расположения точек измерений, а не от их конкретного местоположения в пространстве. Это означает, что корреляционные функции определяются вектором h. Для изотропного случая, когда корреляция не зависит от направления, а только от расстояния, — вектор h переходит в скаляр (расстояние): h = x1 x2.

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

Ковариация (covariance) — статистическая мера корреляции между двумя значениями Z(x) и Z(x + h) в точках, разделенных вектором h:

Для N(h) экспериментальных точек, разделенных вектором h, где m–h — среднее значение для данных, находящихся в началах вектора h:

m+h — среднее значение для данных, находящихся в концах вектора h:

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

Анализ и моделирование пространственной корреляции. Вариография Полувариограмма (semivariogram), или просто вариограмма — вариация разницы значений переменной в двух точках как функция расстояния и направления между ними:

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

Чем ближе значения данных (меньше разница между ними), тем больше значение вариограммы.

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

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

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

• Вариограмма точечно симметрична относительна нуля: (h) = (–h).

• Теоретически в соответствии с определением в нуле вариограмма должна быть равна нулю ((0) = 0), но на практике бывают случаи, когда это условие приходится игнорировать. Они рассмотрены ниже при обсуждении проблем моделирования вариограммы.

• При выполнении условия стационарности второго порядка выполняются следующие соотношения между вариограммой и ковариацией:

где C(0) — априорная ковариация или экспериментальная вариация Var[Z()] функции Z.

В. В. Демьянов, Е. А. Савельева Геостатистика: теория и практика Упражнение 4.1. Вывод выражения для ковариации При условии стационарности среднего значения В соответствии со свойством аддитивности E { A + B} = E { A} + E {B}.

Упражнение 4.2. Связь между ковариацией и вариграммой Показать, что при условии стационарности второго порядка Упражнение 4.3. Симметрия вариограммы Показать, что (h) = (–h).

Существуют статистические моменты, аналогичные вариограмме, но отличающиеся степенью, в которую возводится разница значений пар точек [Goovaerts, 1997].

Мадограмма (madogram) — модуль разницы — позволяет уменьшить влияние больших разбросов значений по сравнению с вариограммой:

Родограмма (rodogram) — квадратный корень — еще более понижает влияние значений с большим разбросом:

Дрейф (drift) — очень важная функция при анализе пространственной корреляции:

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

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

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

Рис. 4.1. Пример поведения характеристик пространственной корреляционной И для ковариации, и для вариограммы существуют стандартизованные варианты: коррелограмма и стандартизованная вариограмма. Эти моменты второго порядка являются более робастными, т. е. более устойчивыми к зашумленным данным и присутствию выбросов (outliers). Они вычисляются по следующим формулам:

коррелограмма (correlogram): ( h) = ;

стандартизованная вариограмма (standardised variogram):

В. В. Демьянов, Е. А. Савельева Геостатистика: теория и практика где –h и +h — стандартные отклонения для точек, находящихся соответственно в начале и в конце вектора h:

Упражнение 4.4. Вариограмма и коррелограмма Определить, при каком предположении верно следующее соотношение:

g ( h) = (0) ( h).

Очень часто в реальных данных наблюдается эффект пропорциональности. Он уже рассматривался в разделе 2.6. Здесь мы рассмотрим проявление эффекта пропорциональности в появлении зависимости между средним m(h) и вариограммой (h). Эффект пропорциональности значительно усложняет понимание результата анализа пространственного корреляционного анализа, внося в вариограмму дополнительную зависимость. Как было описано в Главе 2, обнаружить этот эффект можно с помощью вычисления локальных статистических характеристик — среднего и вариации. Для этого используется метод движущегося окна (см. раздел 2.6). При наличии эффекта пропорциональности предлагается использовать относительные вариограммы, которые нивелируют эффект пропорциональности [Isaaks, Srivastava, 1989].

Можно выделить два вида относительных вариограмм (relative variogram):

• общую относительную вариограмму (general relative variogram):

когда вариограмма просто нормируется на квадрат среднего значения данных, разделенных вектором h;

• парную относительную вариограмму (pairwise relative variogram):

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

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

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

Пространственная ориентация задается вектором h, определяемым длиной (лагом) и направлением. Количество и размер лагов определяются конкретными данными — важно, чтобы было несколько лагов на росте вариограммы и несколько лагов, когда она достигает некоторого уровня, сравнимого со значением априорной вариации, и перестает расти. Если значение вариограммы не перестает расти, это может означать, что для данных не выполнена гипотеза о стационарности второго порядка. Проблема нестационарности более подробно обсуждается в Разделе 4.9.

При подозрении о наличии различий в пространственной структуре в зависимости от направления рассчитываются экспериментальные вариограммы по направлениям (directional variogram). Число направлений обычно определяется количеством данных. Для получения общего представления о наличии анизотропии достаточно двух взаимно перпендикулярных направлений. Для более точного моделирования анизотропной вариограммы удобно иметь 6-8 направлений (см. рис. 4.6). Выбор направления расчета вариограммы чрезвычайно важен для четкого выявления корреляционной структуры.

Исходные данные для анализа обычно произвольным нерегулярным образом распределены по области, поэтому трудно предположить, что удастся набрать достаточное количество пар точек измерений, разделенных точно зафиксированными расстояниями в заданном направлении. Чтобы преодолеть эту проблему, используют допущение по разбросу значения расстояния лага (lag tolerance) и угла раствора вокруг направления (direction tolerance). Допуск расстояния лага Dh (lag tolerance) определяет отклонение расстояния в парах от значения расстояния лага h. При Dh = h/2 все В. В. Демьянов, Е. А. Савельева Геостатистика: теория и практика данные будут учитываться при расчете вариограммы и каждая точка попадет хотя бы в один лаг. Если Dh < h/2, то часть точек может не попасть ни в один лаг из-за ограничения размера допуска расстояния в лаге. Если Dh > h/2, то некоторые данные могут учитываться при расчете значения вариограммы для нескольких лагов. Такое перекрытие лагов бывает полезно при малом количестве данных.

Угол раствора D (direction tolerance) вокруг направления позволяет выявлять узконаправленные анизотропные корреляции. Ширина полосы bw (bandwidth) сужает область поиска на больших расстояниях, ограничивая угол раствора D. Когда угол раствора равен 90°, вариограмма становится обобщенной по всем направлениям (omnidirectional). Такая вариограмма используется, если никакая анизотропия в данных не обнаружена или если анизотропией решено пренебречь, например из-за малого количества данных измерений.

На рис. 4.2 представлена схема параметров для вычисления вариограммы в рамках одного направления : h, Dh, D, bw. Можно представить, что такие сектора (как на рис. 4.2) перемещаются по области данных от одной точки к другой для учета всех пар точек, которые ранжируются по расстоянию между ними и попадают в тот или иной лаг для рассматриваемого направления.

— угол направления вариограммы; D — угол раствора; h — лаг; Dh — разброс лага; bw — (полу)ширина полосы; квадратиками помечены точки измерений, которые используются при вычислении вариограммы Вариограмма, рассчитанная по схеме, изображенной на рис. 4.2, приведена на рис. 4.3а. На графике вариограммы указано число пар для каждого лага. Возрастание вариограммы с расстоянием лага указывает на наличие корреляции между значениями в парах. Скорость роста вариограммы с расГлава Анализ и моделирование пространственной корреляции. Вариография стоянием лага характеризует величины пространственной корреляции. Постоянное значение вариограммы для больших расстояний лага показывает отсутствие корреляции между значениями в парах. Это можно проиллюстрировать диаграммами разброса значений в парах лага (lag scatterplot) для различных лагов, они уже приводились в разделе 2.6. Выберем три лага вариограммы (на рис. 4.3а — 1-й, 2-й и 8-й). Из рис. 4.3б для 1-го лага видно, что значения в парах сгруппированы вдоль диагонали графика — это означает хорошую корреляцию. Однако из-за малого расстояния между точками 1-го лага количество точек для расчета значения вариограммы мало — 52 пары. Для 2-го лага с вдвое большим расстоянием количество пар точек возрастает до 172, что значительно повышает статистическую репрезентативность значения вариограммы. Значения в парах для 2-го лага также имеют высокую корреляцию (рис. 4.3в) в среднем, хотя видно, что разброс значений в парах возрастает. В 8-й лаг вошли пары точек на большом расстоянии (115). На рис. 4.3г видно, что значения в парах для 8-го лага не коррелированны — точки разбросаны по всему графику и не группируются вдоль диагонали, как для 1-го и 2-го лагов.

Рис. 4.3. Вариограмма, количество пар в лагах (а) и диаграммы разброса пар для В. В. Демьянов, Е. А. Савельева Геостатистика: теория и практика Выбор размера отклонений лага и угла направления зависит от количества данных. Если данных много и они распределены плотно, разбросы могут быть небольшими. Важно только следить, чтобы число пар, попавших в каждый сектор, было достаточным. При малом количестве данных допустимо использование перекрывающихся секторов. Они не нарушают общую структуру вариограммы, а делают ее более гладкой, удобной для последующего моделирования. Пример экспериментальных вариограмм, построенных на основе одних и тех же данных, для различных лагов представлен на рис. 4.4. Можно видеть, что с уменьшением значения лага экспериментальная вариограмма становится менее гладкой. Обычно используют равные по длине лаги, что ведет к равномерному их распределению. Однако в особых случаях можно использовать и неравномерно распределенные лаги [Flamm et al., 1994]. Поскольку расчет вариограммы весьма чувствителен к выбору длины лага, их значение является принципиальным для дальнейшего моделирования пространственной корреляции. На практике используют интерактивные программы подбора, такие как описанные в [Pannаtier, 1996;

Kanevski, Maignan, 2004].

Рис. 4.4. Вариограммы, рассчитанные с лагом различной длины:

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

С его помощью можно ограничить разброс направлений пар точек и таким образом выявить узконаправленную корреляцию. На рис. 4.5 приведены вариограммы, рассчитанные для различных значений раствора угла: 45°, 30° и 15°. Можно видеть, что с уменьшением раствора поведение вариограммы становится менее гладким, если направление не полностью соответствует направлению доминирующей пространственной непрерывности.

На рис. 4.5 также видно, что с уменьшением угла раствора число пар, указанное для каждого лага, уменьшается.

Рис. 4.5. Вариограммы по одному направлению для различных углов раствора D:

Для визуализации вариограмм можно использовать графики, как, например, на рис. 4.6, но для изображения и исследования пространственной анизотропии более удобны двумерные изображения. Одним из них является вариограммная роза [Chernov et al., 1998] (рис. 4.7). Роза имеет вид лепестков, представляющих вариограммы по направлениям. Вариограммная роза симметрична относительно центра в силу свойств симметрии моментов второго порядка. По вариограммной розе можно построить изолинии значений вариограммы методом линейной интерполяции на основе триангуляции (см. рис. 4.8). Изолинии вариограмм хорошо визуализируют анизотропную корреляционную структуру.

В. В. Демьянов, Е. А. Савельева Геостатистика: теория и практика Рис. 4.6. Вариограммы по шести направлениям и уровень априорной вариации C(0) Рис. 4.8. Сглаженные изолинии вариограммной розы (светлым контуром отмечен Анализ и моделирование пространственной корреляции. Вариография Другим инструментом, который дает представление о поведении пространственной структуры в целом, является вариограммная поверхность (variogram surface) (рис. 4.9) [Pannatier, 1996]. Для построения вариограммной поверхности вектор h представляется не в полярном виде, как для вариограммной розы (расстояние и направление), а в виде проекций лага на оси координат Dx и Dy. Вариограммная поверхность представляет собой поверхность значений, вычисленных на регулярной сетке в пространстве лагов по формуле вариограммы. При вычислении (наборе пар) вариограммные координаты, естественно, тоже берутся с разбросом (lag tolerance). Вариограммная поверхность позволяет сразу увидеть анизотропию и определить приоритетные направления для построения вариограмм.

Следует отметить, что вариограммная поверхность обладает центральной симметрией относительно точки (0, 0).

Рис. 4.9. Вариограммная поверхность для 137Cs Еще одним инструментом пространственного корреляционного анализа является вариограммное облако (variogram cloud) (рис. 4.10). Это диаграмма разброса квадратов разности значений для всех пар в зависимости от расстояния между точками в паре. Такая диаграмма помогает распознать пары с большим значением квадрата разности значений, поскольку они вносят существенный вклад в значение экспериментальной вариограммы. Присутствие пар, дающих необоснованно большие значения для малых лагов, помогает выявлять крайние экстремальные значения — выбросы (outliers). Вариограммное облако также помогает определить оптимальный лаг для вычисления вариограммы. Вариограммное облако может быть построено для любого направления и раствора угла.

Следует отметить, что вариограммное облако также обладает центральной симметрией относительно точки (0, 0).

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

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

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

Анализ и моделирование пространственной корреляции. Вариография Рис. 4.11. Влияние экстремального значения на вариограмму:

а — вариограммы с учетом крайне высокого значения; б — диаграммы разброса пар для 7-го лага вместе с крайним значением; в — местоположение пар точек с крайним значением; г — диаграммы разброса пар для 6-го лага вместе с крайним значением; д — вариограммы без учета крайнего значения В. В. Демьянов, Е. А. Савельева Геостатистика: теория и практика Упражнение 4.5. Для построения анизотропной модели вариограмм производится расчет вариограммы по нескольким направлениям (углам ).

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

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

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

С другой стороны, система уравнений кригинга имеет единственное решение при несингулярности матрицы системы, что эквивалентно положительной определенности ковариации, что, в свою очередь, эквивалентно отрицательной определенности вариограммы [Armstrong, 1984;

Christakos, 1984]:

где xi — конечное число точек в пространстве (xi: i = 1, 2, 3,..., m); bi — действительные числа (bi, bj: i, j = 1, 2, 3,..., m).

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

Ниже приведены наиболее известные типы моделей вариограмм, удовлетворяющие требованию отрицательной определенности [Goovaerts, 1997].

Модель наггет:

Анализ и моделирование пространственной корреляции. Вариография Константа c0 = C(0) носит название наггет (nugget), что означает самородок. Это понятие было заимствовано из золотодобычи и означает некоррелированный случайный характер. Наличие у данных вариограммы только типа наггет означает отсутствие пространственной корреляции. Данные в этом случае распределены абсолютно случайно (pure nugget). Отсутствие корреляции в данных может иметь следующие причины: мелкомасштабная вариабельность (меньше, чем расстояние между измерениями), ошибки измерений, ошибки в определении местоположений точек.

Сферическая модель:

где a — действительный радиус корреляции (range), на рис. 4.12 a = 40. Радиус корреляции означает, что данные, находящиеся на расстояния a и ближе, коррелированны, а находящиеся на расстоянии больше a — не коррелированны.

Для сферической модели (a) = C(0) = c0 + c — плато (sill). Эта модель ведет себя линейно вблизи нуля.

Экспоненциальная модель:

где а — эффективный радиус корреляции (range), на рис. 4.13 a = 40. На этом расстоянии значение вариограммы достигает 95% плато. Данная модель достигает плато асимптотически.

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

где а — эффективный радиус корреляции (range), на рис. 4.14 a = 40. На этом расстоянии значение вариограммы достигает 95% плато.

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

Степенная модель:

Анализ и моделирование пространственной корреляции. Вариография где — степень. Эта модель (рис. 4.15) при = 1 иногда выделяется в отдельный тип и называется линейной.

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

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

Периодическая модель (hole effect):

где а — период периодической структуры, эквивалентный радиусу корреляции.

Эта модель (рис. 4.16) используется для периодических структур. Периодическая модель работает только в одном направлении.

В. В. Демьянов, Е. А. Савельева Геостатистика: теория и практика Затухающая периодическая модель (dampened Hole effect model):

Эта модель (рис. 4.17) представляет собой произведение экспоненциальной модели ковариации и периодической функции.

Затухающая периодическая структура встречается чаще, чем чисто периодическая.

Рис. 4.17. Затухающая периодическая модель Кубическая модель (Cubic mode):

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

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

Пентасферическая модель (Pentaspherical model):

Эта модель является вариантом сферической модели более высокого порядка. Она также ограничена значением плато c за пределами радиуса корреляции a, на рис. 4.19 a = 50. В отличие от классической сферической модели данная обладает большим градиентом на участке роста.

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

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

Упражнение 4.6. Перечислить стационарные и нестационарные модели вариограмм.

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

Графическое изображение наиболее часто использующихся в геостатистике типов моделей вариограмм приведены на рис. 4.20. Изображенные модели имеют одинаковые параметры: c0 = 0, c = 1, a = 10.

Рис. 4.20. Основные типы стационарных моделей вариограмм при a = 10, c = 1, c0 = Во многих случаях используются линейные комбинации моделей различных типов. Их плато и радиусы могут быть различными. Так, наличие суммы моделей с разными радиусами корреляции свидетельствует о присутствии гнездовой структуры (nested structure) данных:

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

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

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

Индексы качества аппроксимации являются функциями разницы между значением экспериментальной вариограммы (hi) для i-го лага h и значением модели вариограммы *(hi, ) для лага hi и набора параметров модели (с0, c, a). Приведем несколько критериев, которые применяются в геостатистике.

Индекс взвешенных наименьших квадратов:

Суммирование проводится по количеству лагов экспериментальной вариограммы k. Если веса W(i) = 1, то этот метод вырождается в обычный метод наименьших квадратов, который предполагает, что невязки между экспериментальной вариограммой и модельными значениями независимы, нормально распределены и имеют одну и ту же вариацию. Это достаточно сильное предположение, оно не всегда соблюдается на практике.

Индекс Кресси [Cressie, 1985]:

где N(hi) — число пар точек, по которым вычислялось значение для лага.

В. В. Демьянов, Е. А. Савельева Геостатистика: теория и практика Веса для индекса Кресси зависят от количества пар в лаге вариограммы, которое определяется при построении экспериментальной вариограммы (см. раздел 4.3).

Индикатор Кресси может иметь модифицированный вариант [Zhang et al., 1995]:

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

Индекс качества подбора (indicative goodness of fit), используемый в программе VARIOWIN [Pannatier, 1996]:

где hmax(k) — максимальная длина лага для k-го направления; N(hi) — число пар точек, по которым вычислялось значение для лага; n — число лагов;

2 — вариация оценки вариограммы; P — число направлений, которые участвуют в подборе параметров модели.

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

Информационный критерий Акайк [Xiaodong et al., 1996]:

где p — число параметров в модели.

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

4.5. Поведение вариограмм на больших расстояниях Вариограммы приведенных выше типов можно классифицировать по поведению на больших расстояниях. При наличии стационарности второго порядка значение вариограммы на бесконечности равно значению ковариации в нуле (ковариации исходных данных) [Barnes, 1991]. Это характерно для вариограмм сферического, экспоненциального и гауссового типов.

Модель вариограммы сферического типа достигает плато на расстоянии a, экспоненциального типа — на расстоянии 3a, гауссова модель достигает 95%-ного значения плато на расстоянии a 3.

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

Упражнение 4.7. Каково соотношение между ковариацией на бесконечности и вариограммой в предположении стационарности второго порядка?

Упражнение 4.8. Каково соотношение между вариограммой на бесконечности и ковариацией в предположении стационарности второго порядка?

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

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

При этом (0) = 0 в некоторой окрестности нуля, а при h > (h) = C(0), т. е.

чистый наггет-эффект соответствует полному отсутствию корреляций.

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

При линейном поведении вблизи нуля ((h) ~ A|h|) вариограмма не дифференцируема в нуле, но остается непрерывной при h = 0. Этот случай представлен линейной моделью.

4.7. Анизотропия вариограмм До сих пор мы рассматривали модель вариограммы для одного направления или для изотропной вариограммы, которая зависит только от расстояния между точками. При изотропии изолинии вариограммы на вариограммной поверхности или вариограммной розе будут иметь форму окружности. Если вариограмма зависит и от ориентации пары точек в пространстве, то можно говорить о наличии анизотропии. Это означает существование структур данных с различными пространственными характеристиками в различных направлениях. В такой ситуации есть два выхода. Первый выход — построить одну модель изотропной вариограммы и при ее использовании всякий раз производить с вектором h преобразования пространства и только после этого подставлять в качестве аргумента величину |h|. Второй выход состоит в полномасштабном моделировании анизотропной структуры вариограммы и использовании ее в вычислениях. Для моделирования сложной анизотропии используют гнездовую структуру. Для выбора подхода к моделированию анизотропии следует проверить, к какому типу она относится.

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

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

Геометрическая (geometric) анизотропия. В этом случае изолинии вариограммы на вариограммной поверхности или вариограммной розе имеют форму эллипса (рис. 4.21). Это означает, что существует положительно определенная матрица B — такая, что g ( h) = g hT Bh — изотропная вариограмма.

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

где amin — длина малой оси эллипса; amax — длина основной оси эллипса; — угол, определяющий направление основной оси эллипса. Вариограммы по четырем направлениям для случая геометрической анизотропии в горизонтальном направлении (90°) представлены на рис. 4.22. Можно заметить, что графики вариограмм по направлениям следуют в последовательности углов направлений.

Риc. 4.21. Геометрическая анизотропия: вариограммная роза В. В. Демьянов, Е. А. Савельева Геостатистика: теория и практика Рис. 4.22. Геометрическая анизотропия: вариограммы по направлениям Негеометрическая анизотропия радиуса (non-geometric range anisotropy).

В этом случае изолинии на вариограммной поверхности или вариограммной розе образуют форму, отличную от эллипса (пример на рис. 4.23).

При моделировании такого рода вариограмм удобнее пользоваться гнездовыми структурами, хотя при желании можно построить и преобразование пространства [Zimmerman, 1993]. Для такого преобразования сначала нужно построить и запомнить функцию зависимости радиуса корреляции от угла R(), потом необходимо построить изотропную модель вариограммы, взяв за основу одно из направлений. Например, это может быть направление, в котором радиус корреляции максимален и равен R. Тогда при вычислении вариограммы достаточно воспользоваться преобразованием пространства:

где — угол, определяющий направление, в котором достигается максимум значения модуля |h'| (иначе говоря, это угол, для которого вектор h имеет проекцию максимальной длины, т. е. коллинеарный вектору h). Пример, иллюстрирующий негеометрическую анизотропию, представлен на рис. 4.24:

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

Анализ и моделирование пространственной корреляции. Вариография Рис. 4.23. Негеометрическая (зонная) анизотропия радиуса: вариограммная роза Рис. 4.24. Негеометрическая (зонная) анизотропия радиуса: вариограммы Другой случай зонной анизотропии — анизотропия плато (sill anisotropy).

Здесь для различных направлений различаются значения плато (рис. 4.25).

Наличие плато у вариограммы означает, что процесс не только удовлетворяет внутренней гипотезе, но и обладает стационарностью второго порядка. Тогда, учитывая поведение вариограммы на больших расстояниях, для любого фиксированного h можно написать В. В. Демьянов, Е. А. Савельева Геостатистика: теория и практика а — вариограммы по направлениям; б — вариограммная роза Очевидно, что такое равенство справедливо для любого направления, где у вариограммы обнаружено плато. А так как мы имеем различные значения плато для различных направлений, существуют хотя бы два направления h и h2 такие, что lim C ( h1 ) lim C ( h2 ), т. е. по крайней мере в одном направлении lim C ( h ) 0. Таким образом, если значение плато меняется в зависимости от направления, то, следовательно, существует хотя бы одно направление, на котором корреляция между значениями не пропадает ни при каких расстояниях.

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

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

Упражнение 4.9. На рис. 4.26 изображены вариограммы по четырем направлениям: 0°, 30°, 60° и 90°. Найти радиус корреляции в направлениях 270°, 240°, 210° и 180°.

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

В. В. Демьянов, Е. А. Савельева Геостатистика: теория и практика Упражнение 4.10. Вариограммы для различных геологических структур Поставьте пространственным (геологическим) образам на рис. 4.27 в соответствие вариограммы на рис. 4.28 и вариограммные розы на рис. 4.29.

Что можно сказать о корреляционных структурах приведенных образов?

а — одиночное тело; б — извилистое русло; в — эоловые дюны; г — множественные тела; д — параллельные русла; е — вкрапления Рис. 4.28. Вариограммы, соответствующие образам на рис. 4. Анализ и моделирование пространственной корреляции. Вариография Рис. 4.29. Вариограммные розы, соответствующие образам на рис. 4. 4.8. Неоднозначность при моделировании пространственных структур при помощи вариограммы На рис. 4.30а,б приведены два пространственных образа, которые представляют собой совершенно различные геологические структуры: параллельные русла и смыкающиеся эоловые структуры дюн. Эти геологические структуры обладают противоположной связностью — русла связаны горизонтально и способны поддерживать течение потока в себе, в то время как смыкающиеся дюны не имеют сквозной связанности и поэтому не могут поддерживать течение потока. Свойство связности крайне важно при моделировании течения в подземных месторождениях нефти, газа и задачах гидрогеологии.

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

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

Анализ и моделирование пространственной корреляции. Вариография Рис. 4.31. Пространственные образы параллельных русел (а), эллиптических объектов (б), прерывистых русел (в) и соответствующие им вариограммные розы (г, д, е) Для решения подобных задач были разработаны методы многоточечной статистики, базирующиеся на многоточечных структурных функциях, которые можно рассматривать как более общий случай модели пространственной корреляции [Caers, 2005]. Более подробно методы многоточечной статистики изложены в главе 10.

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

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

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

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

Влияние линейного тренда на вариограммы представлено на рис. 4.32.

Видно (см. рис. 4.32а), что данные имеют сильный линейный тренд (отрицательную корреляцию), в результате которого вариограмма ведет себя нестационарно на больших расстояниях, демонстрируя постоянный рост (см. рис. 4.32б). Если вычесть из данных компоненту линейного тренда (y = –0,575x – 0,4712), то оставшиеся невязки (см. рис. 4.32в) демонстрируют стационарную периодическую корреляционную структуру (см. рис. 4.32г), которую не отражала вариограмма данных с трендом. Таким образом, предварительное моделирование тренда и его вычитание из данных позволяют найти локальную корреляцию. Однако на практике модель тренда зачастую является более сложной и нелинейной.

Нелинейный пространственный тренд можно увидеть в форме непрерывной крупномасштабной зависимости наблюдаемой величины от направления (рис. 4.33). Тренд в горизонтальном направлении X можно грубо представить при помощи линейной модели, в то время как тренд в направлении Y имеет ярко выраженный нелинейный характер. При наличии тренда в одном из направлений вариограмма обычно превышает уровень априорной вариации (рис. 4.34). Вариограмма в направлении 45° указывает на присутствие тренда. Вариограмма в направлении 0° обладает стационарностью только на расстоянии до 30. Вариограмма в направлении –45° также указывает на присутствие крайних значений либо тренда.

Анализ и моделирование пространственной корреляции. Вариография Рис. 4.32. Удаление тренда из данных и вариограмма:

а — данные с линейным трендом; б — вариограммы данных с трендом;

в — невязки после удаления линейного тренда; г — вариограммы невязок Рис. 4.33. Зависимости значений переменной в горизонтальном X (а) В. В. Демьянов, Е. А. Савельева Геостатистика: теория и практика Рис. 4.34. Вариограммы по направлениям при наличии нелинейного Присутствие пространственного тренда приводит к невозможности использовать гипотезу о стационарности в каком-либо виде (второго порядка и др.). Следовательно, традиционные геостатистические модели кригинга не могут быть использованы напрямую для интерполяции данных.

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

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

1. Кригинг с трендом (или универсальный — universal — кригинг) использует модель тренда как линейную комбинацию набора базисных функций (см. Раздел 5.4). Универсальный кригинг прост в применении, не требует дополнительных настроек параметров, если правильно выбраны базисные функции. Их выбор и представляет наибольшую трудность. Чаще всего используется полиномиальная модель (линейная комбинация полиномов). Но такая жесткая модель не всегда может адекватно описать сложную многомасштабную пространственную структуру тренда.

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

3. Локально меняющееся среднее (locally varying mean) использует в качестве модели тренда локальное среднее значение, которое может быть получено при помощи метода движущегося окна (см. Главу 2).

4. Кригинг невязок с движущимся окном (moving window residual kriging) похож на модель с локально меняющимся средним. Но он вычислительно гораздо более сложен, поскольку предполагает подбор модели тренда и модели вариограммы в каждой локальной окрестности (окне) [Haas, 1990].

5. Внутренняя случайная функция порядка k (intrinsic random function of order k — IRFk) использует моменты более высокого порядка, чем второй, вместе с вариограммой для моделирования трендов [Marcotte, David, 1988].

6. Моделирование нелинейного тренда на разных масштабах при помощи искусственной нейронной сети (ИНС). Кригинг невязок искусственной нейронной сети (neural netwrok residual kriging — NNRK) был предложен в 1996 г. [Kanevski et al., 1996]. На его основе было разработано целое семейство методов с применением различных типов ИНС (более подробно см. раздел 10.2).

4.10. Пример анализа пространственной корреляционной структуры Для иллюстрации приведем пространственный корреляционный анализ данных по загрязнению поверхности западной части Брянской области изотопом 137Cs (рис. 4.37). Загрязнение произошло после аварии на Чернобыльской АЭС. Оно было принесено облаком, которое частично выпало на поверхность в виде сухого осаждения, а где-то было вымыто локальными дождями.

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

а — диаграмма расположения точек; б — полигоны Вороного Рис. 4.38. Вариограммное облако для данных по загрязнению 137Cs Следующее, на что следует обратить внимание, — возможность присутствия тренда и анизотропии. Для прояснения этих вопросов рассмотрим функцию дрейфа. На рис. 4.39 представлен дрейф, рассчитанный для различных направлений (углы указаны в градусах от направления запад-восток против часовой стрелки). По этому рисунку, а также по розе дрейфа (рис. 4.40) видно, что поведение дрейфа — не колебание вокруг нуля и что оно разГлава Анализ и моделирование пространственной корреляции. Вариография личается в направлениях с северо-востока на юго-запад и с северо-запада на юго-восток. Значит, в данных присутствует тренд, но, судя по значениям модуля дрейфа, незначительный. Для таких данных предпочтительно использовать интерполяционные модели, учитывающие тренд. Но в данном случае традиционные геостатистические методы вполне могут дать приемлемый результат (нужно моделировать направления, в которых дрейф несущественен).

Рис. 4.39. Дрейф по направлениям для данных по загрязнению 137Cs Рис. 4.40. Изолинии розы дрейфа для данных по загрязнению 137Cs В. В. Демьянов, Е. А. Савельева Геостатистика: теория и практика Рис. 4.41. Изолинии розы вариограммы (а) и розы общей относительной вариограммы (б) для данных по загрязнению 137Cs Вариограммная роза (рис. 4.41а) и вариограммная поверхность (рис. 4.42) демонстрируют анизотропную пространственную структуру. В данном случае на вариограммной поверхности и вариограммной розе довольно четко прочертился эллипс. Это указывает на геометрическую анизотропию.

Основная ось расположена вдоль направления 60° от горизонтали, а малая — соответственно –30° от горизонтали.

Рис. 4.42. Вариограммная поверхность для данных по загрязнению 137Cs Теперь можно выполнить моделирование вариограммы, т. е. осуществить подбор математической модели, удовлетворяющей всем свойствам вариограммы и позволяющей описать ее для любого лага и направления (напомним, что экспериментальную вариограмму мы оценили лишь для конечного набора лагов и направлений). Из вариограмм по направлениям наиболее стационарными и подходящими для моделирования являются вариограммы в направлениях 0° и 150°. Они отражают устойчивые корреляционные Анализ и моделирование пространственной корреляции. Вариография структуры данных в этих направлениях. Как уже отмечалось, дрейф в этих направлениях несущественен (см. рис. 4.39), что позволит применить модель обычного кригинга.

Рисунки 4.43, 4.44 иллюстрируют процесс подбора параметров теоретической модели вариограммы при помощи интерактивной программы «Геостат Офис» [Kanevski, Maignan, 2004]. Там же приведены значения индексов (см. Раздел 4.4). Окончательные параметры теоретической модели вариограммы приведены в табл. 4.1. Она имеет гнездовую структуру, состоящую из наггет-модели и двух сферических анизотропных моделей и хорошо аппроксимирует точки экспериментальных вариограмм (рис. 4.45). Для анизотропной модели с геометрической анизотропией был получен радиус корреляции около 60 км вдоль направления 50° (СВ-ЮЗ), что соответствует примерно половине размера области. Зонная анизотропия моделируется второй сферической структурой в направлении 165° (СЗЗ-ЮВВ).

Рисунок 4.46 показывает хорошее соответствие анизотропной структуры экспериментальной и модельной вариограмм.

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

Дальнейшее использование промоделированных вариограмм — решение систем линейных уравнений, что является тривиальным при современном развитии численных методов и вычислительной техники. Основные результаты анализа и их отображение получены с помощью программного обеспечения «Геостат Офис» [Kanevski, Maignan, 1996].

Литература Armstrong M. Common Problems Seen in Variograms // Mathematical Geology. — 1984. — Vol. 16, N 3. — Р. 305—313.

Barnes R. J. The Variogram Sill and the Sample Variance // Mathematical Geology. — 1991. — Vol. 23, N 4. — Р. 673—678.

Caers J. Petroleum Geostatistics / SPE. — [S. l.], 2005. — 98 p.

Chernov S., Demyanov V., Kanevski M., Savelieva E. VarRose — a Way of Variogram Analysis. — Moscow, 1998. — 27 p. — (Препринт / ИБРАЭ;

IBRAE-98-03).

Christakos G. On the problem of permissible covariance and semivariogram models // Water Resources Research. — 1984. — Vol. 20, N 2. — Р. 251—265.

Clark I. Practical Geostatistics. — London; New York: Elsevier Applied Science Publ., 1984.

Cressie N. Fitting models by weighted least squares // Mathematical Geology. — 1985. — Vol. 17, N 5. — Р. 563—586.

David M. Handbook of Applied Advanced Geostatistical Ore Reserve Estimation. — Amsterdam B.V.: Elsevier Applied Science Publ., 1988.

Flamm C., Kanevsky M., Savelieva E. Non-regular variography and multimethod mapping to determination of origin of heavy metals // International Association for Mathematical geology Annual Conference: Papers and Extended Abstracts. — [S. l.], 1994.

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

Oxford Univ. Press, 1997.

Haas T. C. Kriging and automated variogram modeling within a moving window // Atmospheric Environment. — 1990. — Vol. 24A. — Р. 1759—1769.

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

Kanevsky M., Arutyunyan R., Bolshov L. et al. Artificial neural networks and spatial estimations of Chernobyl fallout // Geoinformatics. — 1996. — Vol. 7, N 1—2. — Р. 5—11.

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

Marcotte D., David M. Trend surface analysis as a special case of IRF-k kriging // Mathematical Geology. — 1988. — Vol. 20, N 7. — Р. 821—824.

Pannatier Y. VARIOWIN Software for Spatial Data Analysis. — New York:

Springer-Verl., 1996.

Xiaodong J., Olea R. A., Yu Y.-S. Semivariogram modeling by weighted least squares // Computers and Geosciences. — 1996. — Vol. 22, N 4. — Р. 387—397.

Zhang X. F., Van Eijkeren J. C. H., Heemink A. W. On the weighted least-square method for fitting a semivariogram model // Computers and Geosciences. — 1995. — Vol. 21, N 4. — Р. 605—608.

Zimmerman D. L. Another Look at Anisotropy in Geostatistics // Mathematical Geology. 1993. — Vol. 25, N 4.

Глава Геостатистические интерполяции для одной переменной Данная глава посвящена семейству моделей кригинга для анализа одной пространственной переменной. В Разделе 5.1 формулируются основные постулаты кригинга. Основные типы кригинга (простой и обычный) подробно описаны в Разделах 5.2 и 5.3. В Разделах 5.4 и 5.5 рассмотрены некоторые другие типы кригинга — универсальный, логнормальный. Раздел 5.6 посвящен дополнительным аспектам теории кригинга, в частности некоторым свойствам весовых коэффициентов и вариации кригинга.

Кригинг — базовая интерполяционная модель геостатистики. Он является основой всех методов, связанных с геостатистикой, — интерполяции, вероятностного картирования, стохастического моделирования. Термин «кригинг» служит для обозначения семейства алгоритмов линейной пространственной регрессии. Он происходит от фамилии инженера Д. Крига, который первым применил интерполятор на основе модели пространственной корреляции данных для анализа золотых месторождений Южной Африки [Krige, 1951]. Л. С. Гандин независимо от Д. Крига применил аналогичный метод для объективного анализа метеополей [Гандин, Каган, 1976].

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

5.1. Основные постулаты кригинга Рассмотрим проблему оценивания значения непрерывной переменной Z в произвольной точке x, принадлежащей области пространства S. Исходная информация о переменной представлена в виде набора {z(xi), i = 1, …, n} из n измерений, сделанных в точках x1, x2, …, xn пространства.

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

где i(x) — весовые коэффициенты, относящиеся к данным z(xi). Значения z(xi) интерпретируются как реализации случайной переменной Z(xi).

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

Тип оценивателя зависит от модели случайной функции Z(x). Ее всегда можно разложить на две компоненты — детерминистический тренд m(x) и случайную невязку R(x):

Компонента невязки R(x) моделируется как стационарная случайная функция с нулевым математическим ожиданием mR(x) и ковариацией CR(h):

Математическое ожидание пространственной переменной Z в точке x, таким образом, будет равно значению тренда:

Далее мы рассмотрим разновидности кригинга для моделирования одной переменной, которые определяются предположением о виде тренда.

Все методы семейства кригинга используют одну и ту же целевую функцию для минимизации, а именно вариацию ошибки оценки E ( x ) при дополнительном условии несмещенности оценки, иными словами, вариация минимизируется при ограничении Изначально все кригинги рассматривались как глобальные оцениватели, т. е. для оценки значения в точке x0 из области S предполагалось использовать все имеющиеся измерения {z(xi), i = 1, …, n}. Тогда предположение, например о постоянстве среднего распространяется на всю (возможно, достаточно большую) область, что, вообще говоря, редко встречается в природе. Чтобы не делать такого сильного предположения, на практике обычно используют локальную оценку на основе n(x) ближайших к точке оценивания данных. Можно сказать, что используемые при оценке данные выбираются из некоторой окрестности W(x) точки оценивания x. Размер и форма этой окрестности зависят от исходных данных: предлагается использовать зону, ориентированную в соответствии с эллипсом корреляции, но возможно и меньшего или большего размера. Уменьшение окрестности позволяет получать более вариабельную (менее сглаженную) оценку.

5.2. Простой кригинг Простой кригинг (simple kriging — SK) работает в предположении о стационарности второго порядка случайной переменной Z(x) (см. Раздел 2.6).

Кроме того, предполагается, что детерминистическая компонента m(x) в (5.2) постоянна и известна на всей области исследования S:

Знание среднего значения m дает возможность сделать простое преобразование путем вычета постоянного тренда и далее строить линейный оцениватель для случайной функции Y(x) на всей области S автоматически получая несмещенность оценки (сохранение глобального среднего). Так как E{Y ( x )} = 0, x S, то В. В. Демьянов, Е. А. Савельева Геостатистика: теория и практика Окончательная оценка простого кригинга из (5.5) и (5.6) имеет вид Теперь рассмотрим вариацию ошибки 2 для оценки функции Y(x).

Так как функция Z(x) удовлетворяет стационарности второго порядка, этому условию удовлетворяет и функция Y(x). Тогда все ковариации и вариации, входящие в (5.8), существуют. Чтобы получить вариацию оценки, подставим в первое и второе слагаемые суммы (5.8) формулу оценки (5.6):

Вариация неизвестной случайной переменной Y(x) также существует и связана с априорной вариацией исходных данных 2 :

В итоге получим значение вариации ошибки оценки переменной Y как сумму (5.9), (5.10) и (5.11):

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

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

Оценка функции Z(x) получается подстановкой полученных весовых коэффициентов в формулу (5.7).

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

• Оценка простого кригинга является точной. Это означает, что если координата оцениваемой точки x0 совпадает с какой-то координатой из исходного набора данных (x0 = xi, i = 1,..., n), то полученная оценка будет также совпадать с исходным значением: Z*(x0) = Z(xk). Это легко доказать, пользуясь единственностью решения системы уравнений простого кригинга.

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

• Оценка простого кригинга является сглаженной по сравнению с распределением исходных данных. Как видно из (5.14), вариация оценки простого кригинга 2 меньше значения вариации исходных данных 2.

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

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

5.3. Обычный кригинг Обычный кригинг (ordinary kriging — OK) отличается от простого кригинга тем, что не предполагает знание среднего значения. В обычном кригинге среднее значение считается постоянным, но оно неизвестно.

Кроме того, обычный кригинг при использовании локальной оценки не требует постоянства среднего по всей зоне оценивания; предполагается, что среднее постоянно только в окрестности точки оценивания W(x).

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

Оценка обычного кригинга строится, как линейная комбинация исходных данных:

Рассмотрим условие несмещенности (5.4) в случае неизвестного среднего:

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

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

Решение этой задачи осуществляется с использованием минимизации лагранжиана L(x), куда помимо вариации (5.8) включается условие (5.16) с весовым коэффициентом µ(x) (множителем Лагранжа):

где Cij — ковариации случайных переменных:

Для минимизации лагранжиана L(x) необходимо его продифференцировать по всем весам i(x) и коэффициенту µ(x), а потом приравнять эти производные нулю. В результате получается линейная система из n(x) + 1 уравнений с n(x) + 1 неизвестными — система уравнений обычного кригинга:

Система уравнений (5.17) аналогично с системой уравнений простого кригинга (5.13) имеет единственное решение при положительной определенности функции ковариации С и отсутствии пространственно совпадающих или очень близких точек.

Для вычисления оценки найденные веса i(x) подставляются в линейную комбинацию (5.15). Вариация обычного кригинга вычисляется из формулы (5.12) с использованием первой части системы (5.17):

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

В. В. Демьянов, Е. А. Савельева Геостатистика: теория и практика Вариация обычного кригинга (5.18) также может быть переписана в терминах вариограммы:

Все свойства, описанные в разделе 5.2 для простого кригинга, относятся в той же мере и к обычному кригингу. Hо сравнение формул вариаций простого (5.14) и обычного (5.18), (5.19) кригингов показывает, что платой за неизвестное значение среднего является увеличение вариации, что ведет к росту неопределенности оценки.

Для того чтобы понять, как влияет на оценку и вариацию кригинга отсутствие знания среднего, рассмотрим искусственный пример. Пусть известны одно значение Z(x1) функции Z(x) и модель пространственной корреляции, заданная ковариацией C(h) или вариограммой (h). Используя эти данные, построим оценку функции в точке x0. Если известно среднее (пусть для простоты это будет нуль), модель простого кригинга строится следующим образом:

1. Уравнение простого кригинга:

C11 = C10, где C11 — значение априорной вариации; C10 — значение ковариационной функции для вектора, разделяющего точки x1 и x0.

2. Весовой коэффициент кригинга 3. Определяется оценка:

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

4. Вариация кригинга также определяется пространственной корреляцией.

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

Рис. 5.1. Зависимость оценки и вариации простого кригинга от расстояния при Теперь рассмотрим обычный кригинг, когда среднее неизвестно.

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

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

3. Вариация кригинга 2 = 2 g10 зависит от расположения точек, растет по мере удаления, а потом выходит на удвоенное значение плато.

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

Оценка кригинга для одномерного случая, изображенная на рис. 5.2, демонстрирует точное воспроизведение данных простым и обычным кригингом (по построению оценки). Оценки простого и обычного кригинга весьма близки в приведенном одномерном примере при одинаковых параметрах модели (радиус корреляции r = 50, плато c = 5, наггет c0 = 0). Увеличение наггета до c0 = 2 c соответствующим уменьшением плато до c = 3 ведет к росту случайной компоненты в оценке. Это приводит к более высокой вариации оценки в окрестностях точек измерений вместе с более гладкой оценкой вне этих окрестностей. Таким образом, оценка вне точек измерений стремится к среднему значению, в то время как сами точки измерения воспроизводятся точно в виде «выколотых». Это означает, что данные в модели с высоким наггетом предполагаются зашумленными и менее репрезентативными, т. е. не характерными и поэтому оказывающими слабое влияние вне своей непосредственной окрестности. Вариация соответствующей оценки кригинга (рис. 5.3) демонстрирует, что значения обычного кригинга выше, чем значения простого кригинга. Вариация кригинга с высоким наггетом значительно выше вне точек измерений, чем для модели с нулевым наггетом.

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

Геостатистические интерполяции для одной переменной Рис. 5.2. Оценка простого и обычного кригинга в одномерном случае (для сферической модели вариограммы с радиусом корреляции r = 50 и различных значений наггета с0 = 0 и 2, соответственно плато c равно 5 и 3) Рис. 5.3. Вариация оценки простого и обычного кригинга в одномерном случае Упражнение 5.1. Несмещенность оценки кригинга Показать несмещенность оценки обычного кригинга при условии стационарности переменной Z.

Упражнение 5.2. Оценка кригинга и нестационарность вариограммы Одним из условий кригинга является стационарность (внутренняя гипотеза). Что происходит с оценками кригинга, если нет стационарности среднего значения?

Упражнение 5.3. Точное воспроизведение данных кригингом Кригинг является точным оценивателем — воспроизводит значения исходных данных. Показать, что для заданного значения Z(x0) оценка кригинга Z*(x0) = Z(x0).

Упражнение 5.4. Вариация оценки кригингом Показать, что для заданного значения Z(x0) в отсутствие ошибки измерения вариация оценки кригинга SK(x0) = 0.

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

А. Что больше — вариация оценки кригинга или вариации исходных данных?

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