Какой кригинг выбрать: простой, ординарный, с трендом, блочный, индикаторный

Мы создаем софт для горно-геологических служб калийных рудников, и после первой статьи про кригинг на чистом NumPy самый частый вопрос звучал одинаково: «Хорошо, а какой именно кригинг брать?» Вопрос правильный: под словом «кригинг» живёт целое семейство методов, и выбор между ними влияет на результат сильнее, чем тонкая настройка вариограммы. В плагине Isoliner их пять - простой, ординарный, с полиномиальным трендом, блочный и индикаторный, - и каждый существует не для галочки, а под конкретный класс геологических задач. Под катом - разбор всех пяти на живых примерах калийного месторождения: когда какой, какие параметры за что отвечают, как это выглядит в коде и на каких ошибках спотыкаются чаще всего. Общая механика: что меняется, а что нет Напомню каркас из первой статьи: оценка в точке - взвешенная сумма соседних скважин, веса - решение системы уравнений с ковариациями из вариограммы: s = np.linalg.solve(A, r) # A - ковариации соседей между собой, w = s[:na] # r - ковариации соседей с точкой оценки est = float(np.dot(w, vra)) Все пять видов кригинга решают эту систему. Различаются они тем, что считается неизвестным про среднее поле и что именно оценивается - точка, блок или вероятность. Это и есть ключ к выбору. Простой кригинг (SK): среднее известно Простой кригинг предполагает, что среднее значение поля вам известно заранее и постоянно по площади. Тогда система решается без дополнительных условий, а недобор веса соседей компенсируется этим средним: if ktype == 0: # simple kriging e += (1.0 - w.sum()) * skmean # добор до известного среднего Поведение характерное: вблизи скважин оценка следует данным, а по мере удаления плавно стягивается к заданному skmean . Это одновременно сила и ловушка. Сила - минимальная дисперсия ошибки, когда среднее действительно известно: плотно разбуренный участок со стабильной геологией. Ловушка - «среднее на глаз»: задали 25, а в этом крыле залежи реальное среднее 21, и весь недоразбуренный край уехал на четыре процента вверх. Практическое правило: SK берём, когда среднее подтверждено статистикой по представительной выборке этого же домена. Сомневаетесь - не берите. Простой кригинг: слева среднее задано по данным, справа завышено на 7. Скважины (точки) одни и те же, восточнее пунктира площадь не разбурена. На правом гриде весь недоразбуренный восток «всплыл» к ложному среднему - на семь единиц. Оба грида посчитаны ядром Isoliner (kb2d), синтетика: 42 скважины, сферическая вариограмма. Ординарный кригинг (OK): рабочая лошадь Ординарный кригинг признаёт, что среднего мы не знаем, и оценивает его локально - в каждой окрестности своё. Технически это одно дополнительное уравнение с множителем Лагранжа, требующее, чтобы веса суммировались в единицу: if ktype == 1: # ordinary kriging A[na, :na] = vg.maxcov # условие несмещённости: A[:na, na] = vg.maxcov # сумма весов = 1 A[na, na] = 0.0 r[na] = vg.maxcov Вдали от скважин OK тянется не к глобальному среднему, а к среднему ближайшей окрестности - для месторождений с меняющейся от участка к участку геологией это ровно то, что нужно. Цена - дисперсия ошибки чуть выше, чем у SK с верным средним, но эта осторожность окупается: OK не притворяется, что знает то, чего не знает. Если не знаете, с чего начать, - начните с OK. В Isoliner он и стоит по умолчанию. Кригинг с трендом: когда залежь наклонена У кровли пласта почти всегда есть региональный уклон: моноклиналь, крыло складки. Стационарности - постоянного среднего - нет в принципе, оно закономерно меняется по площади. Кормить такое OK можно, но окрестности будут постоянно «догонять» тренд, а вариограмма получится искажённой: в неё вмешается детерминированная составляющая. Решение - регрессия-кригинг: снять полиномиальный тренд методом наименьших квадратов, кригинговать остатки, вернуть тренд к результату: trend = PolyTrend.fit(xd, yd, vrd, degree) # МНК-полином 1-й или 2-й степени vrd = trend.residuals(xd, yd, vrd) # кригингуем остатки # ... кригинг ... grid += trend.surface(grid_x, grid_y) # тренд возвращается к оценке Isoliner при снятии тренда печатает, какую долю дисперсии он забрал, - если 60-80 процентов, тренд был реальным; если 5 - он вам не нужен. Два предостережения из практики. Первое: вариограмму задавайте по остаткам , а не по исходным данным - плагин напоминает об этом в журнале, потому что наступают все. Второе: полином 2-й степени за пределами облака скважин - парабола, и она улетает; не экстраполируйте грид с квадратичным трендом далеко за данные. Поле с региональным уклоном: слева ординарный кригинг напрямую, справа регрессия-кригинг (плоскость МНК снята, остатки прокригингованы, плоскость возвращена). Восточнее пунктира скважин нет: слева ординарный кригинг глохнет в локальном среднем последних окрестностей, справа снятая плоскость продолжает уклон закономерно. Блочный кригинг: запасы, а не точки Точечный кригинг отвечает на вопрос «сколько в этой точке». Подсчёт запасов спрашивает другое: «сколько в среднем в этом блоке 100 на 100 метров». Средние по блоку спокойнее точечных значений - выбросы внутри блока гасятся, - и блочный кригинг учитывает это математически: правая часть системы усредняется по дискретизации блока, дисперсия ошибки падает. Параметров два: размер блока (обычно шаг будущей эксплуатационной сетки) и дискретизация - на сколько подточек делить блок при усреднении ковариаций (4×4 хватает почти всегда). Одни и те же скважины (с двумя нарочитыми выбросами): слева точечный кригинг, справа блочный с дискретизацией 4×4. Конусы вокруг выбросов на блочном гриде заметно гашены, средняя стандартная ошибка ниже - усреднение по блоку работает. Практическое следствие, которое удивляет новичков: грид блочного кригинга глаже точечного, и это не потеря информации, а корректный переход к другому носителю. Сравнивать содержание пробы и содержание блока - классическая ошибка носителя (support effect), из-за которой запасы «теряются» на бумаге. Индикаторный кригинг: категории и вероятности Минтип, фация, зона замещения - категориальные признаки, их не проинтерполируешь напрямую. Индикаторный кригинг превращает категорию в индикатор (1 внутри класса, 0 вне), кригингует эти нули-единицы обычным OK и получает вероятность класса в точке . Дальше либо карту вероятности используют как есть (карта риска замещения - готовый продукт для планирования), либо режут порогом на домены. Тонкости: вариограмма строится по индикаторам (у неё свой, обычно более короткий радиус - в наших данных по зонам замещения ~3.7 км против 6-8 у содержания), а значения вне [0, 1], которые кригинг может дать у резких границ, принудительно обрезаются. Индикаторный кригинг: скважины с категорией (красные - замещение, белые - сильвинит) превращаются в карту вероятности P(замещение). Порогом 0.5 из неё получается контур домена, а сама карта - готовый слой риска для планирования. Окрестность поиска: общие параметры всех пяти Какой бы вид ни выбрали, окрестность - радиус поиска, минимум и максимум соседей, анизотропия - влияет на результат не меньше. Три правила, снимающие большинство проблем: Радиус порядка радиуса вариограммы. Меньше среднего расстояния между скважинами - появятся дыры «недостаточно соседей»; в разы больше радиуса корреляции - лишние вычисления без пользы. Максимум соседей 12-16. Больше - система тяжелее, а дальние точки всё равно получают исчезающие веса. Анизотропия окрестности согласована с анизотропией вариограммы: вытянутая залежь - вытянутый эллипс поиска, азимут и отношение осей берутся прямо из вариограммной карты (о ней была первая статья). Полный пример: грид в двадцать строк Ядро kb2d не зависит от QGIS - это чистый NumPy, поэтому всё из статьи запускается обычным питоном. Вот код, которым посчитаны картинки выше, целиком: import numpy as np from grid_isolines.kb2d import Variogram, build_grid rng = np.random.default_rng(11) x = rng.uniform(40, 620, 42) # скважины: запад разбурен, y = rng.uniform(40, 960, 42) # восток - нет v = 25 + 4*np.sin(x/170) + 3*np.cos(y/220) + rng.normal(0, 0.8, 42) vg = Variogram(nugget=2.5, # C0 и структура - в абсолютных structures=[dict(it=1, # единицах; it: 1 сферическая, cc=12.0, # 2 экспоненциальная, 3 гауссова, aa=210.0, # 4 степенная ang=0.0, anis=1.0)]) grid = build_grid(x, y, v, vg, ktype=1, skmean=0.0, # 0 = SK (нужен skmean), 1 = OK ndmin=1, ndmax=14, # окрестность: мин/макс соседей rad2=420.0**2, # радиус поиска, в квадрате nodata=-999.0, xmn=0, ymn=0, cell=1000/72, nx=72, ny=72) Дальше по вкусу: ndisc=4 превращает вызов в блочный кригинг, with_variance=True возвращает вторым гридом стандартную ошибку, индикаторный - это тот же OK по нулям-единицам категории. Модуль берётся из установленного плагина (папка grid_isolines в профиле QGIS) или из клона репозитория - лицензия GPL, пользуйтесь. Шпаргалка Задача Вид кригинга Ключевой параметр Содержание внутри однородного домена, данных много Простой (SK) Обоснованное среднее Универсальный случай, старт любой задачи Ординарный (OK) Радиус и вариограмма Кровля/подошва с региональным уклоном С трендом Степень 1, редко 2 Запасы по сетке блоков Блочный Размер блока Минтип, замещение, любая категория Индикаторный Порог и своя вариограмма Попробовать за пять минут Всё выше воспроизводится без собственных данных: инструмент «Создать пример скважин (демо)» генерирует облако с трендом, анизотропией и категориальным минтипом - гоняйте все пять видов на одном наборе и сравнивайте гриды и кросс-валидацию. В следующей статье отдельно разберём изолинии: от маршевых квадратов до бергштрихов и контурных полигонов, границы которых совпадают с линиями геометрически точно. Ссылки: Плагин: https://plugins.qgis.org/plugins/grid_isolines/ Код: https://github.com/Valery35/qgis-isoliner (GPL) Первая статья цикла: https://habr.com/ru/articles/1055612/ Вопросы и рабочие случаи приветствуются - особенно спорные: «какой кригинг тут был правильным» разбирать интереснее всего. Isoliner развивается на задачах реальных предприятий. Если вашему производству не хватает функции - напишите нам: https://www.informpp.ru/ Только зарегистрированные пользователи могут участвовать в опросе. Войдите , пожалуйста. Сколько видов кригинга вы планируете использовать или протестировать: 0% 1 0 0% 2-3 0 100% 4-5 1 Проголосовал 1 пользователь. Воздержавшихся нет.