Мы создаем софт для горно-геологических служб калийных рудников, и после первой статьи про кригинг на чистом 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), синтетика: 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 хватает почти всегда).
Практическое следствие, которое удивляет новичков: грид блочного кригинга глаже точечного, и это не потеря информации, а корректный переход к другому носителю. Сравнивать содержание пробы и содержание блока - классическая ошибка носителя (support effect), из-за которой запасы «теряются» на бумаге.
Индикаторный кригинг: категории и вероятности
Минтип, фация, зона замещения - категориальные признаки, их не проинтерполируешь напрямую. Индикаторный кригинг превращает категорию в индикатор (1 внутри класса, 0 вне), кригингует эти нули-единицы обычным OK и получает вероятность класса в точке. Дальше либо карту вероятности используют как есть (карта риска замещения - готовый продукт для планирования), либо режут порогом на домены.
Тонкости: вариограмма строится по индикаторам (у неё свой, обычно более короткий радиус - в наших данных по зонам замещения ~3.7 км против 6-8 у содержания), а значения вне [0, 1], которые кригинг может дать у резких границ, принудительно обрезаются.
Окрестность поиска: общие параметры всех пяти
Какой бы вид ни выбрали, окрестность - радиус поиска, минимум и максимум соседей, анизотропия - влияет на результат не меньше. Три правила, снимающие большинство проблем:
Радиус порядка радиуса вариограммы. Меньше среднего расстояния между скважинами - появятся дыры «недостаточно соседей»; в разы больше радиуса корреляции - лишние вычисления без пользы.
Максимум соседей 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://github.com/Valery35/qgis-isoliner (GPL)
Первая статья цикла: https://habr.com/ru/articles/1055612/
Вопросы и рабочие случаи приветствуются - особенно спорные: «какой кригинг тут был правильным» разбирать интереснее всего.
Isoliner развивается на задачах реальных предприятий. Если вашему производству не хватает функции - напишите нам: https://www.informpp.ru/
