Мы создаем софт для горно-геологических служб калийных рудников, и это третья статья цикла про геостатистику в QGIS: в первой был кригинг на чистом NumPy, во второй — выбор между пятью его видами. Сегодня две темы, которые остались за кадром, а в повседневной работе геолога решают едва ли не больше самой интерполяции. Первая — изолинии: не картинка поверх растра, а документ, который идёт в подсчёт запасов и на печатные планы, со строгими требованиями к геометрии. Вторая — всё, что живёт вокруг кригинга и превращает интерполяцию в решения: карта стандартной ошибки, кросс-валидация, внешний дрейф, вероятность превышения кондиции, гидравлика потока по Дарси и гауссова симуляция неопределённости.
Изолинии: почему «просто contour» не подходит
Стандартный путь — gdal_contour или родные изолинии QGIS — даёт линии. Геологу нужны линии и контурные полигоны между ними («пояса»): закрашенные интервалы значений, по которым меряются площади и считаются запасы. И вот тут начинается геометрия: границы полигонов обязаны совпадать с изолиниями точно, без зазоров и самопересечений, а пояса — покрывать сто процентов площади грида. Если линии и полигоны строить двумя независимыми проходами, совпадения не будет никогда: разные алгоритмы, разные допуски.
Isoliner решает это в одну сторону: изолинии трассируются маршевыми квадратами по гриду, а пояса получаются разрезанием плоскости этими же линиями — QgsGeometry.unaryUnion плюс polygonize из GEOS. Полигоны рождаются из линий, поэтому их границы тождественны по построению. Забавная производственная деталь: сначала пояса резались через штатный алгоритм native:splitwithlines, но в QGIS 4 он начал терять кусочки покрытия, и прямые вызовы GEOS оказались и точнее, и быстрее.
Скриншот QGIS: изолинии с поясами на демо-гриде, панель слоёв с линиями и полигонами.
Сглаживание: форма против точности
Сырые маршевые квадраты дают ломаные с изломом на каждой ячейке — технически верно, на печатном плане неприемлемо. В плагине два независимых механизма. Бикубическая денсификация уплотняет линию по гриду — изолиния следует полю точнее, чем исходная ломаная. Скругление Чайкина срезает углы чисто геометрически — форма приятнее, но линия слегка отходит от поля. Оба на чистом NumPy, оба с параметрами, и правило простое: денсификация почти всегда, Чайкин — когда план смотрят люди, а не алгоритмы, и отход на доли ячейки не имеет значения.
Бергштрихи
Депрессионные изолинии без бергштрихов — коротких штрихов вниз по склону — российский геолог на план не примет: без них замкнутая изолиния не говорит, впадина внутри или холм. Плагин ставит штрихи автоматически по градиенту грида, со стороной и шагом в параметрах. Мелочь, которой нет почти нигде, а без неё «почти готовый» план отправляется на доработку.
Скриншот QGIS: замкнутая депрессионная изолиния с бергштрихами.
Где гриду верить: карта стандартной ошибки
Кригинг — единственный из популярных интерполяторов, который вместе с оценкой отдаёт дисперсию ошибки в каждой точке. Isoliner пишет её вторым выходом — гридом стандартной ошибки, и читать его надо так: ошибка минимальна у скважин, растёт в пустотах сети и взлетает за границей разбуренности. Это карта качества разведочной сети: где сгущать скважины, докуда продлевать контур подсчёта, где грид — оценка, а где уже экстраполяция.
Одно важное «но» из первой статьи: стандартная ошибка живёт в масштабе вариограммы. Если силл задан «на глаз», карта ошибки будет систематически врать — и проверяется это кросс-валидацией.
Кросс-валидация: MSDR и QQ
Скользящий контроль leave-one-out: каждая скважина по очереди выбрасывается и предсказывается по остальным. Из остатков считаются два диагноза. MSDR — средний квадрат стандартизованного остатка: при согласованной модели он равен единице, меньше — ошибки переоценены (силл завышен), больше — недооценены. QQ-график стандартизованных остатков против нормальных квантилей показывает то же самое глазами: точки на диагонали — модель согласована, наклон меньше единицы — переоценка ошибок.
Практический цикл подбора выглядит так: экспериментальная вариограмма → первая модель → кросс-валидация → MSDR далёк от единицы, крутим наггет и силл → повторяем. Обычно сходится за два-три оборота, и после этого карте стандартной ошибки можно верить буквально — в тех же единицах, что и данные.
Кригинг с внешним дрейфом: тренд, который не полином
Во второй статье тренд снимался полиномом от координат. Но бывает, что «тренд» — это не гладкая функция x и y, а известное всюду поле другой величины. Мощность пласта закономерно связана с отметкой кровли, содержание — с сейсмическим атрибутом, уточняемая поверхность — с грубой региональной моделью по редкой сетке. Кригинг с внешним дрейфом добавляет к скважинам это знание о форме поля между ними: снимается линейная регрессия на растр-дрейф, кригингуются остатки, регрессия возвращается к оценке. Та же схема регрессия-кригинг, только объясняющая переменная берётся из грида, а не из формулы.
Требований к дрейфу два: он должен быть известен во всех точках будущего грида (иначе там не будет оценки) и связь должна быть реальной — плагин печатает в журнал долю дисперсии, которую забрала регрессия, и если это единицы процентов, дрейф выдуман, возвращайтесь к обычному OK.
Скриншот QGIS: окно «Кригинг с внешним дрейфом», растр соседней поверхности как дрейф.
Вероятность превышения: из оценки и ошибки — в решения
Подсчёт запасов и планирование спрашивают не «сколько», а «какова вероятность, что содержание выше кондиции». Раз кригинг отдаёт в каждой ячейке оценку и стандартную ошибку, ответ считается одной формулой: локальное распределение принимается нормальным, и вероятность превышения порога — это функция нормального распределения от расстояния оценки до порога в единицах ошибки. Где оценка уверенно выше кондиции — вероятность около единицы, где ниже — около нуля, а вдали от скважин, где ошибка велика, карта тянется к 0.5 — и это правильное поведение: там, где данные молчат, уверенности нет ни в какую сторону.
Инструмент работает постобработкой: берёт готовые растры оценки и ошибки, поэтому одинаково обслуживает обычный кригинг и кригинг с дрейфом.
Градиент и направление потока: геометрия фильтрации
Для гидрогеолога поверхность напора из кригинга — это полдела, дальше нужен ответ «куда и насколько круто течёт». Инструмент градиента берёт растр напора и отдаёт три слоя: модуль градиента (безразмерный i = |∇h| — тот самый гидравлический градиент), азимут направления потока (ноль — север, по часовой) и прореженный точечный слой векторов, который сразу оформлен стрелками: маркер повёрнут по полю азимута, размер масштабирован по модулю градиента, картина фильтрации читается без единого клика в символике. Направление считается строго: вода течёт вниз по градиенту, от большего напора к меньшему, а на пологих участках, где направление не определено, азимут остаётся пустым — вместо выдуманной стрелки будет её отсутствие.
Практическая тонкость, о которой забывают: дифференцирование усиливает шум. Даже аккуратный грид кригинга после взятия градиента способен дать пятнистое поле с дёргаными стрелками, поэтому в инструменте есть лёгкое сглаживание напора перед расчётом (в ячейках грида) — включают его не для красоты, а чтобы вернуть полю градиента читаемость. Важно понимать границы: этот инструмент — геометрия потока, не скорость. Он отвечает «куда» и «насколько круто», но не «как быстро».
Закон Дарси: от геометрии к физике
Скорость появляется, когда к градиенту добавляется проницаемость. Удельный расход по Дарси: q = K·|∇h|, метры в сутки — объём воды через единицу площади сечения. Если вместо коэффициента фильтрации задать растр водопроводимости, считается расход через единицу ширины потока Q = T·|∇h| в квадратных метрах в сутки — и поскольку T это K, умноженный на мощность, мощность пласта отдельно не нужна. Истинная скорость движения воды — это q, делённое на эффективную пористость. Пористости в данных обычно нет, инструмент её не запрашивает, а при необходимости деление делается растровым калькулятором.
Откуда берутся растры K и T: свойства известны в точках опытных откачек и наливов, всюду их разносит тот же кригинг. И здесь главная тонкость раздела — логарифм. Коэффициент фильтрации и водопроводимость почти всегда лог-нормальны, их значения разбросаны на порядки, и кригинг по самим значениям исказит поле: пара высокопроницаемых точек перетянет всю окрестность. Правильно кригинговать логарифм. В «2D Kriging» для этого есть преобразование ln: кригуется ln(K), а растр возвращается уже в исходных единицах — тогда в инструменте Дарси ничего дополнительно включать не нужно. Если же растр приготовлен заранее в логарифмах, флаг «Растры K и T заданы как ln» вернёт значения экспонированием. И ещё одно правило из практики: напорные и безнапорные пласты кригуйте раздельно — у них разная физика мощности, смешивать их в одной вариограмме нельзя.
Весь путь проходится на демо: у скважин из генератора «Создать пример скважин» есть поля head, K и T — постройте напор кригингом по head, подайте его в градиент, добавьте прокригингованный в логарифмах K, и цепочка «уровни в скважинах → напор → поток → расход» замкнётся, не покидая плагина и не требуя ни одного собственного файла.
Гауссова симуляция: сколько на самом деле неопределённости
Кригинг даёт одну сглаженную поверхность. Последовательная гауссова симуляция (SGS) отвечает на другой вопрос: строит ансамбль равновероятных реализаций, каждая из которых воспроизводит гистограмму и вариограмму данных, проходит через скважины и потому остаётся шершавой, как настоящая геология. По ансамблю в каждом узле набирается распределение значений: среднее реализаций сходится к кригингу, а разброс показывает неопределённость напрямую — без предположения о нормальности, которое было нужно карте превышения.
Механика: значения переводятся в нормальные баллы, узлы обходятся в случайном порядке, в каждом простой кригинг по соседям и уже симулированным узлам даёт локальное распределение, из которого разыгрывается значение — и тут же становится обуславливающим для следующих. Скважины зафиксированы во всех реализациях.
Код выпуска: MSDR и вероятность превышения
Как и во второй статье, всё считается чистым NumPy без QGIS. Скользящий контроль и главный диагноз модели — десять строк:
import numpy as npfrom grid_isolines.kb2d import Variogram, cross_validate# x, y, v - скважины из примера второй статьиvg = Variogram(nugget=0.3, structures=[dict(it=1, cc=1.8, aa=180.0, ang=0.0, anis=1.0)])est, var = cross_validate(x, y, v, vg, ktype=1, skmean=0.0, ndmin=1, ndmax=14, rad2=420.0**2, nodata=-999.0)ok = (est != -999.0) & (var > 0)z = (v[ok] - est[ok]) / np.sqrt(var[ok]) # стандартизованные остаткиprint("MSDR =", np.mean(z**2)) # цель - около единицы
Крутите наггет и силл, пока MSDR не придёт к единице, — это и есть цикл подбора из текста. Параметры выше уже подобраны им (MSDR 0.99), а с силлом 14.5 «с запасом» тот же код печатает 0.15.
Карта вероятности превышения из готовых гридов оценки и ошибки — одна формула через функцию ошибок:
from math import erf, sqrtgrid, se = build_grid(x, y, v, vg, ktype=1, skmean=0.0, with_variance=True, ndmin=1, ndmax=14, rad2=420.0**2, nodata=-999.0, xmn=0, ymn=0, cell=1000/72, nx=72, ny=72)thr = 26.0P = 0.5 * (1.0 + np.vectorize(erf)((grid - thr) / (se * sqrt(2.0))))
Этими двумя фрагментами посчитаны QQ-пара и светофорная карта из статьи — воспроизводимость полная.
Шпаргалка выпуска
-
Пояса и изолинии обязаны рождаться из одной геометрии — иначе не сойдутся никогда.
-
Денсификация — точность, Чайкин — внешний вид. Включаются независимо.
-
Карта стандартной ошибки — это карта качества сети, а не украшение.
-
MSDR → 1 — главный винт подбора наггета и силла. QQ показывает то же глазами.
-
Внешний дрейф — когда тренд известен как поле, а не как формула. Долю снятой дисперсии смотрим в журнале.
-
Вероятность превышения переводит оценку и ошибку на язык решений: брать, не брать, доразведать.
-
Градиент — геометрия потока. Дифференцирование усиливает шум, лёгкое сглаживание напора возвращает читаемость.
-
Дарси — физика: q = K·i, Q = T·i. K и T лог-нормальны, кригуйте их логарифм (преобразование ln в «2D Kriging»).
-
SGS — когда нужна неопределённость без предположений: ансамбль вместо одной поверхности.
Всё воспроизводится на демо-данных плагина за пять минут: «Создать пример скважин» → кригинг со стандартной ошибкой → кросс-валидация → изолинии с поясами. В следующий раз — самое зрелищное: собственный 3D-просмотр, тела пластов из многоканальных гридов и блочная модель, которая опрашивается кликом.
Ссылки:
-
Код: https://github.com/Valery35/qgis-isoliner (GPL)
Расскажите про свои требования к изолиниям-документам — у маркшейдеров и гидрогеологов они заметно разные, и нам интересно, чего не хватает вашим планам.
Isoliner развивается на задачах реальных предприятий. Если вашему производству не хватает функции — напишите нам: https://www.informpp.ru/главная-страница/предприятиям
ссылка на оригинал статьи https://habr.com/ru/articles/1057232/