Работа с анизотропией при блочном моделировании посредством интерполяторов кригинга, IDW и DVWE

 

В.А.Мальцев

 

ВНИИГЕОСИСТЕМ, 113105, Россия, Москва, Варшавское шоссе 8

E-mail:  vl-maltsev@yandex.ru.

 

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

ABSTRACT. The main aim of the paper is to discuss technique of modeling of the anizotripy when using kriging or IDW estimators. Special study was done on typical mistakes, related to mixing between the anisotropy of the ore and the anisotropy of the sampling grid. Recommendations are carried out and illustrated with case studies. Several improvements were carried out to the IDW estimators, allowing to model the anisotropy more carefully. This chain of improvements leads to a completely new estimator, called DVWE, that is variogram-based, but uses math, similar to IDW. The discussion section is on comparing this estimator with kriging and IDW, both theoretically and via case studies.

Введение.

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

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

Экспериментально-иллюстративная часть выполнена не вынесенными в общий пользовательский интерфейс средствами программы Geostatistical Software Tool. Большая часть приводимых примеров двумерна, но это только для удобства восприятия статьи, при переходе к трехмерным моделям практически ничего не меняется.

Случай, когда интерполятором является кригинг

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

Рассмотрим два последних утверждения.

Со вторым всё совсем просто. Система уравнений, к примеру, «обычного» кригинга выглядит как

  ,

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

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

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

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

       

Рис.2. Справа – пример типового распределения весов проб при кригинге, слева – пространственное взаиморасположение оцениваемого блока, проб в первых окружениях и более удаленных, а также вариантов геометрии поискового эллипсоида.

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

1. Декластеризация. Разведочная сеть практически никогда не бывает равноплотной во всех направлениях. Ситуация, когда расстояние между профилями метров 200, между скважинами в профиле метров 40, а шаг опробования вдоль скважины 2 метра, достаточно типовая. Существует масса методов для «уравнивания» такой неравномерности опробования. Правильный выбор поискового эллипсоида просто является еще одним таким же методом, причем одним из самых эффективных.

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

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

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

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

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

Если при применении кригинга вариограмму всё же моделируют, в том числе с анизотропией и неправильный выбор эллипсоида портит результат только когда он уж очень неправильный, то в случае интерполяции обратными расстояниями (IDW) ситуация радикально меняется. Хотя большинство программ и имеют инструменты для достаточно корректного учета анизотропии на уровне алгоритма, инструменты эти достаточно глубоко запрятаны, моделирование вариограммы не производится, а в результате в большинстве случаев анизотропию пытаются смоделировать исключительно выбором поискового эллипсоида. Это в корне неверный подход. IDW, так же, как и кригинг, является весовым методом, и, точно так же, как и в кригинге, веса проб за вторым окружением становятся малы. Не исчезающее малы, об этом ниже, но достаточно малы, чтобы определенная через эллипсоид анизотропия в очень значительной степени «замазывалась» ближними к центру блока пробами, конфигурация которых остается неизменной при любой геометрии эллипсоида. Тем самым, все вышеизложенные соображения сохраняются и здесь. Посмотрим, как с ними следует бороться.

1. Вариант «традиционный», имеющийся в большей части программ. Анизотропия задается как новая система координат, преобразованием в которую является трехмерный поворот и сжатие (или растяжение) по двум из трех осей. Направления координат этой новой системы являются направлениями анизотропии, коэффициенты растяжения координатных осей – аффинными модулями. Некоторые программы вместо двух модулей требуют задать три радиуса, но это одно и то же – аффинные модули являются просто отношениями пар этих радиусов. Суть здесь в том, что преобразование сводит изменчивость к изотропной. А далее в метод IDW просто вместо исходного расстояния между каждой парой точек подается расстояние между ними в этой новой системе координат. Разумеется, для понимания направлений и коэффициентов анизотропии придется воспользоваться вариограммой или любой другой структурной функцией, но достаточно на визуальном уровне, без функционального её описания, как это требуется для кригинга.

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

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

2. Второй вариант, лишенный последнего недостатка – это IDW с ограничением радиусов влияния. Для этого нужно добавить к описанию анизотропии всего лишь один параметр – предельный радиус влияния по главному направлению (тому, по которому сжатие/растяжение координатной оси отсутствует). Математически же данная модификация осуществляется тривиально: вместо подстановки расстояния , где  - направления новых осей, а - модули растяжения/сжатия с первой оси на вторую и третью, используем подстановку

.

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

Рис.3. Различные варианты учета анизотропии при интерполяции методом обратных квадратов расстояний. 1 – анизотропия учтена преобразованием расстояний через направления и модули растяжения (модуль с X на Y 1.5, модуль с X на Z 4.0), а эллипсоид взят шаром радиуса 60, 2 – анизотропия моделируется геометрией эллипсоида, выровненного по осям анизотропии с полуосями 60:40:15), 3 – анизотропия определена как в 1, но эллипсоид чисто декластеризующий, выровнен по скважинам и уплощен по длине композитов (полуоси 100:100:15).

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

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

Дальнейшая модификация метода обратных расстояний.

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

Рис.4. Пример трёххструктурной вариограммы (четыре направления – темные линии, осредненная – серая) с разной симметрией структур.

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

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

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

превращается в интерполятор

     ,

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

.

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

Посмотрим пример сопоставления на том же массиве результатов работы кригинга и построенного алгоритма (назовём его DVWE, то есть интерполятор с прямым взвешиванием учетом вариограммы, direct variogram weighting estimator).

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

Легко видеть, что алгоритм DVWE не только гораздо интереснее работает с анизотропией, но и вообще по свойствам чрезвычайно близок к обычному кригингу. Меньший локальный контраст, чем в вариантах IDW, главным образом связан с тем, что IDW не обрабатывает эффекта самородков, который здесь достаточно высок (0.3).

Зачем нужен интерполяционный алгоритм DVWE

Может возникнуть вполне естественный вопрос – зачем нужен еще один кригингообразный интерполятор, лишённый при этом главных достоинств кригинга, а именно экранного эффекта и почти доказуемой оптимальности?

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

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

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

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

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

Краткие выводы.

1. Моделирование анизотропии при блочном моделировании не является тривиальной задачей и имеет свои особенности при применении различных типов интерполяторов.

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

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

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

4. Дальнейшим развитием интерполятора IDW, гораздо более корректно работающим с анизотропией, является предлагаемый в данной статье интерполятор DVWE.

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

О литературе.

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