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

Мы создаем софт для горно-геологических служб калийных рудников, и после первой статьи про кригинг на чистом 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 берём, когда среднее подтверждено статистикой по представительной выборке этого же домена. Сомневаетесь — не берите.

Оба грида посчитаны ядром Isoliner (kb2d), синтетика: 45 скважин, сферическая вариограмма.

Ординарный кригинг (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 хватает почти всегда).

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

Индикаторный кригинг: категории и вероятности

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

Тонкости: вариограмма строится по индикаторам (у неё свой, обычно более короткий радиус — в наших данных по зонам замещения ~3.7 км против 6-8 у содержания), а значения вне [0, 1], которые кригинг может дать у резких границ, принудительно обрезаются.

Окрестность поиска: общие параметры всех пяти

Какой бы вид ни выбрали, окрестность — радиус поиска, минимум и максимум соседей, анизотропия — влияет на результат не меньше. Три правила, снимающие большинство проблем:

• Радиус порядка радиуса вариограммы. Меньше среднего расстояния между скважинами — появятся дыры «недостаточно соседей»; в разы больше радиуса корреляции — лишние вычисления без пользы.

• Максимум соседей 12-16. Больше — система тяжелее, а дальние точки всё равно получают исчезающие веса.

• Анизотропия окрестности согласована с анизотропией вариограммы: вытянутая залежь — вытянутый эллипс поиска, азимут и отношение осей берутся прямо из вариограммной карты (о ней была первая статья).

Полный пример: грид в двадцать строк

Ядро kb2d не зависит от QGIS — это чистый NumPy, поэтому всё из статьи запускается обычным питоном. Вот код, которым посчитаны картинки выше, целиком:

import numpy as npfrom grid_isolines.kb2d import Variogram, build_gridrng = 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, пользуйтесь.

Шпаргалка

Попробовать за пять минут

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

Ссылки:

• Плагин: 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/