Геостатистика в QGIS без SAGA: кригинг на чистом NumPy

Мы создаем софт для горно-геологических служб калийных рудников. Наши геологи и маркшейдеры каждый день превращают тысячи скважинных проб в карты: отметки кровли пласта, содержания KCl, мощности, газоопасность. Классический инструмент для этого — кригинг, и в QGIS он формально есть: SAGA, GRASS, Smart-Map, связки со SciPy. На практике же каждый вариант чем-то не устраивал, и год назад я начал писать свой плагин. Сейчас Isoliner — это 24 инструмента в официальном репозитории plugins.qgis.org: кригинг четырёх видов, вариограммный анализ, кросс-валидация с отчётами, изолинии с контурными полигонами, геологические разрезы и собственный 3D-просмотр. Вычислительное ядро — чистый NumPy, ни одной внешней зависимости.

Под катом — зачем понадобился ещё один кригинг, как выглядит система кригинга в двадцати строках NumPy, что такое вариограмма на пальцах и почему абсолютные единицы силла — главные грабли всех, кто впервые открывает такие инструменты.

Зачем ещё один кригинг

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

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

Кросс-валидация — не роскошь. Модель вариограммы, подобранная на глаз, обязана быть проверена: leave-one-out, распределение остатков, MSDR. В большинстве плагинов этого этапа просто нет, а без него кригинг превращается в генератор красивых, но ничем не подтверждённых картинок.

Скользкие зависимости. Официальный репозиторий QGIS не позволяет объявлять pip-зависимости, а просить геолога открыть OSGeo4W Shell и набрать pip install — верный способ потерять пользователя. Значит, ядро должно работать на том, что уже есть в любой сборке QGIS: NumPy и GDAL.

QGIS 4. Плагин обязан жить и в QGIS 3.16, и в QGIS 4 на Qt6 — а это отдельная сага о переехавших enum’ах и изменившихся API.

Русский язык. Наши пользователи работают на русском; интерфейс, справка каждого инструмента и 80-страничное руководство должны быть двуязычными.

Ядро: GSLIB KB2D на NumPy

За основу вычислительного ядра взят KB2D — двумерный кригинг из GSLIB, классической стэнфордской библиотеки геостатистики, описанной у Deutsch & Journel. Алгоритм переписан на векторный NumPy, но математика — каноническая.

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

Сердце — сборка и решение системы для каждой точки грида:

# A — ковариации между соседями, r — ковариации соседей с точкой оценкиfor i in range(na):    for j in range(i, na):        c = vg.cova2(xa[i] - xa[j], ya[i] - ya[j])        A[i, j] = c        A[j, i] = c    r[i] = rhs[i]if ktype == 1:                    # ординарный кригинг: несмещённость    A[na, :na] = vg.maxcov        # добавляем множитель Лагранжа    A[:na, na] = vg.maxcov    A[na, na] = 0.0    r[na] = vg.maxcovs = np.linalg.solve(A, r)w = s[:na]                        # веса соседейest = float(np.dot(w, vra))       # оценка

Плюс окрестность поиска (радиус, максимум соседей, эллипс анизотропии), четыре модели вариограммы (сферическая, экспоненциальная, гауссова, степенная), варианты SK/OK, блочный и индикаторный кригинг — но принципиально всё сводится к этой системе. Никакого SciPy: np.linalg.solve покрывает всё, а для устойчивости гауссовой модели принудительно удерживается минимальный наггет — иначе матрица становится почти вырожденной и по гриду расползаются осцилляции.

Скорость на практике: грид 500×500 по нескольким тысячам скважин считается за десятки секунд на обычном офисном компьютере — для ежедневной работы геолога этого достаточно, а предсказуемость чистого NumPy дороже процентов производительности.

Вариограмма на пальцах

Всё, что кригинг «знает» о геологии, он знает из вариограммы, поэтому половина инструментов Isoliner — про неё.

Полувариограмма γ(h) описывает, насколько статистически связаны значения в двух точках в зависимости от расстояния h между ними: для каждой пары точек берётся половина квадрата разности значений, эти величины усредняются по интервалам расстояния — и получается кривая. Чем меньше γ, тем теснее связь.

У кривой три характеристики:

• Наггет C0 — значение, к которому кривая стремится при h → 0. Теоретически там должен быть ноль (точка сравнивается сама с собой), но на практике остаётся ступенька: погрешность измерений плюс изменчивость на масштабе мельче шага сети. Наггет управляет характером поверхности: C0 = 0 — кригинг точный интерполятор, обязан пройти через каждую скважину, и одиночный выброс превращается в конус («бычий глаз»); C0 > 0 — оценка вблизи скважин подтягивается к локальному среднему, поверхность сглаживается.

• Силл — уровень полки, на которую выходит кривая: насколько в среднем различаются далёкие друг от друга скважины. Практически равен дисперсии данных.

• Радиус a — расстояние выхода на полку; дальше него точки статистически не связаны.

И вот здесь живут главные грабли, на которые наступает каждый первый пользователь любого геостатистического софта: наггет и силл задаются в абсолютных единицах дисперсии данных, а не в долях от единицы. Значение силла=1 по умолчанию — условность, которую почти всегда нужно менять. Живой пример по содержанию KCl: среднее ≈ 25 %, дисперсия ≈ 47.6 %², значит, суммарный силл ≈ 47.6, а наггет при желаемой доле 0.35 — это 17, а не «0.35». На сам грид абсолютный масштаб не влияет (важно лишь отношение C0 : C), но карта стандартной ошибки и MSDR в кросс-валидации будут в реальном масштабе только при силле ≈ дисперсии. Isoliner печатает дисперсию данных в журнал при каждом запуске — это и есть ориентир.

Отдельный инструмент строит вариограммную карту: γ по всем направлениям сразу, полярная диаграмма. Если залежь вытянута — корреляция вдоль простирания дальше, чем поперёк, на карте виден эллипс, и его азимут с отношением осей идут прямо в параметры анизотропии кригинга.

Кросс-валидация: доверяй, но проверяй

Подобрали модель — проверяем: leave-one-out по всем скважинам, каждая по очереди выбрасывается и предсказывается по остальным. Инструмент выдаёт HTML-отчёт: гистограмма остатков, QQ-график стандартизованных остатков, MSDR (средний квадрат стандартизованного остатка при удачной модели ≈ 1), и карту худших остатков — точки, где модель промахивается сильнее всего, часто это либо ошибки в данных, либо граница геологического домена.

Рабочий цикл получается коротким: экспериментальная вариограмма → подбор модели → кросс-валидация → если MSDR далёк от единицы или QQ загибается, крутим наггет и радиус → кригинг. Параметры между инструментами переносятся профилями обработки, чтобы не перенабирать по десять полей.

Изолинии, которые сходятся с полигонами

Финальный продукт геолога — план изолиний с закрашенными поясами. Isoliner строит изолинии маршевыми квадратами по гриду, сглаживает их (бикубическая интерполяция плюс скругление Чайкина — оба на NumPy), а пояса получает не повторной трассировкой, а разрезанием плоскости самими изолиниями: QgsGeometry.unaryUnion + polygonize из GEOS. Границы полигонов и линии совпадают тождественно, покрытие — 100 % площади. Забавная деталь: сначала пояса резались через native:splitwithlines, но в QGIS 4 этот алгоритм начал терять кусочки покрытия, и прямые вызовы GEOS оказались не только точнее, но и быстрее.

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

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

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

Что дальше

За кадром остались геологические разрезы (полный конвейер от линии на карте до чертежа с угловыми точками и 3D-забором) и собственный 3D-просмотрщик с телами пластов — про него будет отдельная история: как затащить pyqtgraph и PyOpenGL внутрь плагина и пройти сканер безопасности plugins.qgis.org, который прогоняет Bandit по 8672 файлам и блокирует релиз за один shell=True в чужом коде.

Ссылки:

• Плагин: https://plugins.qgis.org/plugins/grid_isolines/

• Код: https://github.com/Valery35/qgis-isoliner (GPL)

• Руководство (RU/EN, PDF) — в комплекте плагина, открывается из окна «О плагине»

Вопросы и критика геостатистиков приветствуются — особенно про то, чего не хватает вам в повседневной работе с кригингом в QGIS.