У геологов свой майнкрафт: как построить то, что не знаешь, по тому, что знаешь


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


    Посчитаем, сколько там нефти


    Чтобы посчитать, сколько нефти находится в месторождении, как она распределена в пространстве, и как будет идти процесс её добычи, умные люди строят разного рода модели месторождения. Исходной информацией для моделирования являются, в первую очередь, различного рода исследования стволов скважин. Также привлекают множество другой дополнительной информации, которая косвенным или прямым образом характеризует породу или насыщающие её жидкости (нефть, воду или газ). Ну а чтобы жизнь нефтяного инженера не была похожей на сказку, исходные данные имеют различную точность, различную физическую природу, распределены в пространстве неравномерно и характеризуют разный объём породы.


    Существуют разные модели, позволяющие охарактеризовать месторождение. Например, если мы по скважинам узнаем среднюю толщину продуктивного пласта $h$, по изъятому керну посчитаем среднюю пористость $\phi$ и нефтенасыщенность $s$, «прикинем» площадь месторождения $S$, то формула $h\cdot\phi\cdot s\cdot S$ является простейшей (можно назвать нуль-мерной) моделью запаса нефти. (Если слова пористость и нефтенасыщенность вызывают у вас дискомфорт, приглашаем прочитать одну из предыдущих статей нашего блога)


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


    Пусть месторождение описывается трехмерной областью $V$, а $u$ — это точка из множества $V$. Тогда пусть $\phi\left(u\right)$ будет значением пористости в точке $u$, а $s\left(u\right)$ — значением нефтенасыщенности. Тогда формула $\intop_{V}s\left(u\right)\cdot\phi\left(u\right)du$ — будет уже трехмерной моделью запаса нефти. Задачей инженера является построение трехмерных полей $\phi\left(u\right)$, $s\left(u\right)$ (а также проницаемости, насыщенности газа, воды и других величин) в объёме месторождения $V$. А перед этим конечно нужно ещё оконтурить сам геометрический объект $V$. И всё это можно назвать задачами интерполяции, подразумевая, что на основе точечных измерений тех или иных пространственных величин в скважинах нужно построить прогноз распределения их значений между скважинами (на самом деле не только между ними, но также и в неразбуренной области, в этом смысле это и задачи экстраполяции тоже).


    В этой статье мы попробуем кратко и популярно поговорить о методах интерполяции в приложении к геологическому моделированию. С полным осознанием того, что каждая формула по Стивену Хокингу уменьшает количество читателей в два раза, а к данному моменту было приведено уже две формулы, нам всё же не удастся добраться до конца, не написав ещё несколько формул. Мы сделаем это в максимально понятной форме — не для того вы уже прочли столько слов, чтобы уйти со страницы, ничего не поняв. К тому же цель нашего блога — обозначить проблематику и показать актуальные задачи нефтянки для специалистов в других областях, чтобы зародить диалог и, мы надеемся, взаимовыгодное сотрудничество. Поэтому призываем, физически оставаясь на диване, мысленно с него приподняться и двинуться в дебри нефтяной инженерной науки.


    Осложнять нам путь будет множество факторов. Попробуем некоторые из них перечислить; это, собственно, и будут те вызовы, с которыми имеют дело нефтяные инженеры и разработчики софта для нефтяного и геологического моделирования.


    Косвенный характер исходных данных


    Поскольку информация о породе на глубине в несколько километров стоит очень дорого, в ход идёт всё, что может прямо или косвенно охарактеризовать наш продуктивный пласт. Напомним, нефтяников интересуют, прежде всего, свойства пустотного пространства, а не сама порода. В первую очередь, это упомянутые пористость и проницаемость.


    Однако непосредственно пористость и проницаемость мы можем лицезреть лишь в редких образцах породы, извлекаемых при бурении разведочных скважин с отбором керна (как он выглядит, мы показывали в этой статье). Объём таких исследований очень невелик и явно недостаточен, чтоб описать месторождение. К тому же на стол лаборатории попадают только цилиндры размером в несколько сантиметров, высверленные в выборочных местах из колонки керна, который не рассыпался при извлечении и лежит в ящиках достаточно большими и прочными кусками.


    Поэтому сами данные керна в интерполяции и прогнозе межскважинного пространства не участвуют, а участвуют данные, полученные косвенным образом в виде корреляций через посредство геофизических исследований, которые уже проводят практически на всех скважинах. В скважину после бурения спускают гирлянду геофизических приборов, которые записывают сигналы различной физической природы, которые позже сопоставляются со значениями керна. Полученные закономерности используются уже во всем фонде скважин, чтобы рассчитать интересующие нас параметры породы. Этим занимается геофизическая наука. Если вам кажется, что в такой косвенности нет никакой проблемы, то представьте себе, что в травмпункте кто-то пытливый посчитал корреляцию между плотностью костной ткани и ростом пациентов с переломами, а потом применил эту зависимость ко всему населению.


    Итак, косвенный характер исходных данных — первый осложняющий фактор, который надо иметь в виду и осознавать его последствия.


    Носитель или объём осреднения


    Второй особенностью является объём осреднения, который неявно стоит за любым значением, таким как пористость или проницаемость и т.д. Если сделать мысленный эксперимент и рассмотреть точку пласта, в ней может находиться либо камень, либо жидкость. Самой пористости в точке нет, а возникает она в результате объёмного осреднения, когда рассматривается определённого размера и формы окрестность точки, в ней рассчитывается объём пор и делится на объём этой области. Область называется математически носителем или — в англоязычной литературе — словом «support». Тоже самое, только с дополнительными сложностями, относится и к проницаемости, поскольку надо понимать, проницаемость в каком направлении имеется в виду. Безотносительно носителя говорить об объёмных величинах не имеет никакого смысла.


    Южноафриканский учёный Денни Криге, обрабатывая статистические данные запасов полезных ископаемых, экспериментально установил зависимость вариации величин от объёма осреднения (Krige, D.G. 1951. A statistical approach to some basic mine valuation problems on the Witwatersrand. Journal of the Chemical, Metallurgical and Mining Society of South Africa, December 1951. pp. 119–139.). Что, в общем-то, очевидно, чем на большие куски мы мысленно распилим породу, тем меньше стоит ожидать изменчивость величин, которые усредняются по их объёму. Значит, пористость, измеренная в кусочках керна размером 5×5 см, и пористость, оценённая по геофизическим исследованиям скважин, охватывающим пласт породы глубиной до 1 м, и пористость в ячейках цифровой геологической модели размером 50 м×50 м×1 м — это всё разные пористости, и у них разные статистические закономерности, связанные друг с другом правилами перемасштабирования (ещё одно иностранное слово — upscaling).


    Хватит воды


    Давайте теперь уже напишем несколько формул. Пусть мы говорим об интерполяции некоторой пространственной величины $z\left(u\right)$, где $u$ — это точка в пространстве или на плоскости. И пусть в некоторых точках нам известны значения $z_{i}=z\left(u_{i}\right)$. Тогда функция многих аргументов


    $\bar{z}\left(u\right)=F\left(u,z_{1},..,z_{n},u_{1},...,u_{n}\right)$


    будет интерполятором. Естественно требовать, чтобы выполнялись тождества


    $z_{k}\equiv F\left(u_{k},z_{1},..,z_{n},u_{1},...,u_{n}\right), (1)$


    то есть интерполяция в точку с уже известным значением должна давать это самое значение.


    Чаще всего строятся линейные интерполяторы вида


    $\bar{z}=\sum_{i}z_{i}\nu_{i},$


    где $\nu_{i}$ — некоторые веса, которые зависят от точек $u$ и $u_{i}$. И откроем программистам и айтишникам великую тайну, почему математики так любят линейные модели. А любят они их не потому, что они лучше описывают мир, а потому что это единственное, что зачастую они могут посчитать.


    Давайте напишем самый простой интерполятор:


    $\bar{z}\left(u\right)=\sum_{i}z_{i}\frac{1}{\left|u-u_{i}\right|^{m}}\cdot\left(\sum_{i}\frac{1}{\left|u-u_{i}\right|^{m}}\right)^{-1},(2)$


    где $m$ — некоторый показатель. Если $m=1$ то это метод обратных расстояний (см. рис. 1). А если $m=2$, то это метод обратных квадратов расстояний. Метод обратных расстояний будет давать острые пики в точках данных, метод обратных квадратов расстояний будет давать гладкие перегибы в точках данных. Очевидно в эту формулу нельзя просто так подставить в правую часть $u=u_{k}$, потому что возникнет деление на ноль. Но зато, если постепенно приближать $u$ к $u_{k}$, вес при соответствующем $z_{k}$ будет стремиться к 1, в то время как веса при всех остальных $z$ будет стремиться к нулю, поэтому $\bar{z}\left(u_{k}\right)$ принудительно приравнивается к $z_{k}$.


    Ниже приведена простейшая Python программа для интерполяции методом обратных расстояний.


    Inverse distance interpolation
    import numpy as np
    import matplotlib.pyplot as pl
    
    # num of data
    N = 5
    
    np.random.seed(0)
    
    # source data
    z = np.random.rand(N)
    u = np.random.rand(N)
    
    x = np.linspace(0, 1, 100)
    y = np.zeros_like(x)
    
    # norm weights
    w = np.zeros_like(x)
    
    # power
    m = 2
    
    # interpolation
    for i in range(N):
        y += z[i] * 1 / np.abs(u[i] - x) ** m
        w += 1 / np.abs(u[i] - x) ** m
    
    # normalization
    y /= w
    
    # add source data
    x = np.concatenate((x, u))
    y = np.concatenate((y, z))
    order = x.argsort()
    x = x[order]
    y = y[order]
    
    # draw graph
    pl.figure()
    
    pl.scatter(u, z)
    pl.plot(x, y)
    
    pl.show()
    
    pl.close()


    Рисунок 1. Интерполяция методом обратных расстояний


    Самосогласованность


    Проведём мысленный эксперимент. Предположим, мы воспользовались определённым алгоритмом интерполяции, чтобы предсказать глубину пласта в какой-либо точке на основе ранее пробуренных скважин в окрестности. Предположим затем, что мы пробурили в данном месте скважину, и прогноз был полностью оправдан. Какая радость и для геолога-инженера, и для разработчика софта, и для прикладного математика, что их прогноз оправдался, верно? Это убедит нас, что предположения приняты верно и используемая математическая модель неплоха и должна использоваться дальше. Давайте теперь включим пробуренную скважину в список исходных данных и применим ту же математическую модель, только теперь данных у нас больше на 1 штуку. Логично ожидать, что прогноз не должен никак измениться, ведь предыдущая скважина только подтвердила его. А вот на те, этот тест проваливают примитивные алгоритмы пространственной интерполяции, такие как описанный ранее метод обратных расстояний (см. рис. 2). Поэтому такие алгоритмы из арсенала интерполяционных методов сразу отправляются на свалку.



    Рисунок 2. Несамосогласованность метода обратных расстояний


    Декластеризация


    Существует такой хрестоматийный кейс, когда горный инженер, обрабатывая статистические данные некоторого рудного месторождения, провел осреднение содержания полезной руды в образцах, умножил на объём месторождения и получил чрезвычайно привлекательные цифры. Только он не учёл, что преимущественные места отбора были сосредоточены как раз в самых лакомых участках месторождения, то есть данные его несколько тенденциозны. Такой ошибки не должны совершать и нефтяники, поскольку первые скважины очевидно бурятся в самых продуктивных местах. Декластеризация — имя этому действию, смысл которого в том, чтобы придать данным определённый вес в соответствии с их расположением: данные, находящиеся в непосредственной близости друг от друга, значат меньше, чем данные в отдельных точках. Насколько меньше и как это считать, вытекает из статистической модели пространственных величин.


    Ещё немного математики


    Взглянем на формулу (2) немного под другим ракурсом:


    $\bar{z}(u)=\sum_{i} z_i\cdot f_i(u), (3)$


    то есть результат можно понимать, как линейную комбинацию некоторых функций координат, каждая из которых равняется единице в соответствующей точке данных и нулю во всех остальных. Так как функции построены таким специальным образом, нам не пришлось решать никаких систем линейных уравнений, поскольку $i$-ая функция никак не влияет на $k$-ую точку. Но эта простота и делает алгоритм несамосогласованным.
    Построим другой, более сложный интерполятор:


    $\bar{z}(u)=\sum_{i}\lambda_i\cdot c\left(∥u-u_i∥\right), (4)$


    где функция от расстояния $c(h)$ как правило монотонно-убывающая и стремящаяся к нулю на бесконечности, а $\lambda_i$ — весовые коэффициенты. То есть модель (4) — это линейная комбинация одинаковых функций «шапочек», центры которых расположены в точках данных. Осталось правильно рассчитать веса $\lambda_i$. Что бы сделать это, надо добиться выполнения тождеств (1):


    $z_k=\sum_{i}\lambda_i\cdot c\left(∥u_k-u_i∥\right), (5)$


    что порождает систему линейных уравнений на $\lambda$. В матричном виде мы можем записать её в виде $Z=C\cdot\lambda$, где $Z$ — столбец значений $z_k$, $\lambda$ — столбец искомых весов, $C$ — матрица взаимных значений функции «шапочки» $c(h)$ между точками данных. Математик тут сразу запишет, что $\lambda=C^{-1}\cdot Z$, и подставит это в (6), в итоге получит что


    $\bar{z}(u)=\sum_{i}z_i\cdot g_i\left(u\right), (6)$


    где


    $g_i(u)=\sum_{k}C^{-1}_{i, k}\cdot c\left(∥u_k-u∥\right). (7)$


    Формула (6) является одним из вариантов интерполяции методом кригинга (имени того самого Дени Криге). А формула (4) является другим вариантом записи того же кригинга, только в слегка развёрнутом виде, именуемым в англоязычной литературе как dual kriging. Сравните формулы (6) и (3), вроде тоже самое, только функции $g_i$ построены иначе. Если подставить вместо $u$ точку $u_k$, то в формуле (7) опытный глаз заметит, что матрица, обратная к матрице $C$, умножается на один из её столбцов, из чего следует что $g_i(u_i)=1$ и $g_i(u_k)=0$ для любых $i≠k$ (так же как у функций $f_i$).


    В любой литературе по численному моделированию геологии вы неизбежно увидите подобные формулы и встретитесь с интерполяцией методом кригинга (см. рис. 3). А всё потому, что эта интерполяция самосогласована и выполняет декластеризацию, о которой мы говорили выше.


    Предлагаю убедиться, что кригинг — это всего несколько программных строк на языке Python:


    Kriging interpolation
    import numpy as np
    import matplotlib.pyplot as pl
    
    # num of data
    N = 5
    
    np.random.seed(0)
    
    # source data
    z = np.random.rand(N) - 0.5
    u = np.random.rand(N)
    
    x = np.linspace(0, 1, 100)
    y = np.zeros_like(x)
    
    # covariance function
    def c(h):
        return np.exp(-np.abs(h ** 2 * 20.))
    
    # covariance matrix
    C = np.zeros((N, N))
    
    for i in range(N):
        C[i, :] = c(u - u[i])
    
    # dual kriging weights
    lamda = np.linalg.solve(C, z)
    
    # interpolation
    for i in range(N):
        y += lamda[i] * c(u[i] - x)
    
    # add source data
    x = np.concatenate((x, u))
    y = np.concatenate((y, z))
    order = x.argsort()
    x = x[order]
    y = y[order]
    
    # draw graph
    pl.figure()
    
    pl.scatter(u, z)
    pl.plot(x, y)
    
    pl.show()
    
    pl.close()


    Рисунок 3. Интерполяция методом кригинга


    Сглаживающий эффект


    Ничего принципиально более совершенного, чем интерполяция методом кригинга для геологии, до сих пор не придумано. Именно формулы (5) и (6) применяются в дорогущих пакетах геологического моделирования, когда строятся пространственные поля пористости и проницаемости для оценки запасов месторождения.


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


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


    Интерполяция на основе имеющихся данных даёт нам наиболее ожидаемое значение в каждой точке пласта. Для подсчёта запасов нефти это как раз то, что нужно. Однако для построения гидродинамической модели это нам принципиально не подходит. О том, что придумано взамен, пойдёт речь дальше.


    Детерминистика и стохастика


    Поднимем немного температуру на нашем математическом градуснике и импортируем модуль теории вероятности и линейной алгебры.


    Libraries
    from theory import probability
    from numpy import linalg

    У кого данные библиотеки не установлены, наверное, придётся ограничиться неполным пониманием текста.


    Чтобы преодолеть сглаживающий эффект метода кригинга, было придумано стохастическое моделирование. Стохастический подход рассматривает распределение наших пространственных величин в пласте как стохастическое (случайное) поле: в каждой точке пласта u мы рассматриваем $z(u)$ как случайную величину, но все эти случайные величины связаны друг с другом (не независимы друг от друга). И задача стоит синтезировать один или несколько вероятных представителей случайного поля при том, что нам известны некоторые значения $z_i=z(u_i)$.


    Заметим здесь важное логическое отличие. Интерполяция даёт нам оценку или наиболее ожидаемое значение в каждой точке пласта. Стохастическая же модель претендует на синтез таких распределений свойств пласта в межскважинном пространстве, которые могли бы быть похожими на реальность (с точки зрения принятых предположений и гипотез, конечно).


    Само по себе утверждение, что что-то под землёй является случайным — не более чем математическая абстракция. Но чтобы что-то начать считать, нужно эту модель ещё конкретизировать. Математик однажды сказал следующее: пусть в каждой точке пространства $z(u)$ — гауссова случайная величина (имеющая нормальное распределение вероятности), а в совокупности любой набор случайных гауссовых величин $z(u_i)$, взятых в произвольном наборе точек $u_i$, является мультигауссовым, то есть подчиняется специальной многомерной функции распределения. То есть $z(u)$ — это мультигауссово случайное поле, реализацию которого мы хотим синтезировать.


    Для любых двух мультигауссовых случайных величин их ковариация (ковариация по определению $Cov\left(z_1, z_2\right)=E\left(\left(z_1-E\left(z_1\right)\right)⋅\left(z_2-E\left(z_2\right)\right)\right)$, её можно образно сравнить с углом между двумя векторами в пространстве) исчерпывающим образом описывает их взаимозависимость. А взаимозависимость набора случайных величин задаётся ковариационной матрицей — ковариацией каждой случайной величины с каждой. А для случайного поля, в котором в любой точке u соответствует своя случайная величина $z(u)$, нужно построить ковариационную функцию.


    Математик дополнительно предположил, что случайное поле геологического параметра является стационарным, то есть ковариационная функция является функцией только расстояния между точками: $Cov\left(z\left(u\right),z\left(v\right)\right)=c\left(∥u-v∥\right)$. Да-да, ковариационная функция $c(h)$ — это та самая функция «шапочка» в интерполяции кригинга в формуле (4), но здесь более полно раскрывается её смысл: $c(h)$ количественно описывает возможность прогнозировать значение параметра на расстоянии $h$ по его значению в точке. Например, если $c(500)/c(0) = 0.5$, то на расстоянии $500$ метров наша уверенность в прогнозе составляет $50$ у.е., а если $c(1000) = 0$, то знание о значении в данной точке абсолютно бесполезно для прогноза того, что будет через километр.


    Стационарная гауссова случайная модель приобрела большую популярность не потому, что она лучше описывает мир, а потому, что это практически единственное, что мы можем посчитать. Дело в том, что гауссовы случайные величины можно умножать на числа и складывать друг с другом, при этом будут получаться другие гауссовы случайные величины. И ещё про них доказано, что если ковариация гауссовых случайных величин равна нулю, то они независимы. Сказанное позволяет любой мультигауссовый набор случайных величин $Z$ выразить в виде линейной комбинации независимых гауссовых случайных величин $\zeta$ с нулевым математическим ожиданием и единичной дисперсией. В матричном виде запишется


    $Z=A⋅ζ.$


    Чтобы найти матрицу $A$, нужно выразить получившуюся ковариационную матрицу (попарную ковариацию каждого элемента $Z$ с каждым). В матричном виде:


    $C=E\left(Z\cdot Z^T\right)=E\left(A\cdot\zeta\cdot\zeta^T\cdot A^T\right)=A\cdot E\left(\zeta\cdot\zeta^T\right)\cdot A^T=A\cdot A^T,$


    здесь мы воспользовались тем, что случайные величины $\zeta$ независимы, и, значит, ковариационная матрица $E\left(\zeta\cdot\zeta^T\right)$ является единичной. Последняя формула даёт нам ключ к стохастическому моделированию. Для набора случайных величин $Z$, которые мы хотим синтезировать, нужно составить ковариационную матрицу $C$, найти любое её разложение вида $C=A\cdot A^T$ (таких много) и тем самым выразить набор взаимосвязанных случайных величин $Z$ через набор независимых случайных величин $\zeta$, а это то, что математику и нужно. Последние легко можно синтезировать с помощью генератора псевдослучайных величин. Это называется стохастическим гауссовым моделированием.


    Примечательно, что стохастическое моделирование в самом своём принципе подразумевает возможность моделировать вообще без исходных данных (а если быть точным, то без самих скважинных данных, но с предположениями о неоднородности моделируемого геофизического параметра пласта — среднее, дисперсия и вариограмма). Ниже приведём простейшую программу стохастического гауссова моделирования, не учитывающую исходные данные (см. рис. 4).


    Unconditional stochastic gaussian modeling
    import numpy as np
    import matplotlib.pyplot as pl
    
    np.random.seed(0)
    
    # source data
    N = 100
    x = np.linspace(0, 1, 100)
    
    # covariance function
    def c(h):
        return np.exp(-np.abs(h ** 2 * 250))
    
    # covariance matrix
    C = np.zeros((N, N))
    
    for i in range(N):
        C[i, :] = c(x - x[i])
    
    # eigen decomposition
    w, v = np.linalg.eig(C)
    
    A = v @ np.diag(w ** 0.5)
    
    # you can check, that C == A @ A.T
    
    # independent normal values
    zeta = np.random.randn(N)
    
    # dependent multinormal values
    Z = A @ zeta
    
    # draw graph
    pl.figure()
    
    pl.plot(x, Z)
    
    pl.show()
    
    pl.close()


    Рисунок 4. Несколько реализаций стохастического гауссового процесса без учёта данных


    Для учёта исходных данных в стохастическом моделировании нужно строить постериорную ковариацию или, иными словами, нужно построить условную функцию распределения (при условии, что нам известны значения в некоторых точках). Это уже немного сложно объяснить в двух словах, но для любопытных мы приведём пример исходного кода и результат на рисунке 5.


    Conditional stochastic gaussian simulation
    import numpy as np
    import matplotlib.pyplot as pl
    
    np.random.seed(3)
    
    # source data
    M = 5
    
    # coordinates of source data
    u = np.random.rand(M)
    
    # source data
    z = np.random.randn(M)
    
    # Modeling mesh
    N = 100
    x = np.linspace(0, 1, N)
    
    # covariance function
    def c(h):
    return np.exp(−np.abs(h ∗∗ 2 ∗ 250))
    
    # covariance matrix mesh−mesh
    Cyy = np.zeros ((N, N))
    for i in range (N):
        Cyy[ i , : ] = c(x − x[i])
    
    # covariance matrix mesh−data
    Cyz = np.zeros ((N, M))
    
    # covariance matrix data−data
    Czz = np.zeros ((M, M))
    for j in range (M):
        Cyz [:, j] = c(x − u[j])
        Czz [:, j] = c(u − u[j])
    
    # posterior covariance
    Cpost = Cyy − Cyz @ np.linalg.inv(Czz) @ Cyz.T
    
    # lets find the posterior mean, i.e. Kriging interpolation
    lamda = np.linalg.solve (Czz, z)
    y = np.zeros_like(x)
    
    # interpolation
    for i in range (M):
        y += lamda[i] ∗ c(u[i] − x)
    
    # eigen decomposition
    w, v = np.linalg.eig(Cpost)
    A = v @ np.diag (w ∗∗ 0.5)
    
    # you can check, that Cpost == A@A.T
    
    # draw graph
    pl.figure()
    for k in range (5):
        # independent normal values
        zeta = np.random.randn(N)
        # dependent multinormal values
        Z = A @ zeta
        pl.plot(x, Z + y, color=[(5 − k) / 5] ∗ 3)
        pl.plot(x, Z + y, color=[(5 − k) / 5] ∗ 3, label=’Stochastic realizations’)
    
    pl.plot(x, y, ’. ’, color=’blue’, alpha=0.4, label=’Expectation(Kriging)’)
    pl.scatter(u, z, color=’red ’, label=’Source data’)
    pl.legend()
    pl.show()
    pl.close()


    Рисунок 5. Несколько реализаций стохастического гауссова процесса с учётом данных


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


    Если сделать $100500$ стохастических реализаций геологической модели, посчитать среднее значение между ними в каждой точке, то окажется, что мы получим очень близкую картину к результату интерполяции методом кригинга (6). Таким образом, данный метод даёт нам математическое ожидание или самое вероятное значение в каждой отдельно взятой точке пласта, а стохастическое моделирование даёт нам равновероятные случайные картины распределения геологического параметра во всём объёме пласта сразу (каждая из которых, между прочим, совсем не похожа на интерполяцию методом кригинга).


    Попытаемся привести ещё одно рассуждение, поясняющее разницу между интерполяцией и моделированием. Предположим, мы пытаемся сделать прогноз температуры на июль. Очевидно в средней полосе России должно быть тепло и в целом солнечно. Однако каждый знает, что погода у нас редко бывает долго стабильная, и теплые дни неизменно сменяются прохладными, а солнце — дождями и грозами. Например, в среднем после пяти солнечных дней должно быть два пасмурных. Синтез такого правдоподобного погодного графика — это моделирование. А утверждение, что в среднем будет что-то около 22 градусов — интерполяция. Аналогичная интерполяция на март даст вам -5 С, а моделирование — постоянные прыжки туда-сюда через ноль. Очевидно, что тот, кто будет пытаться на основании этого делать предсказания на вредное влияние климата на асфальт, сделает диаметрально противоположные выводы!


    А как происходит моделирование в пространстве?


    Приведённые выше примеры являются одномерными и тривиальными. В случае геологического моделирования пространственные объекты (2D и 3D), где все теоретические выкладки имеют тот же математический смысл (что и для 1D случая), а, с точки зрения программирования, разница только в размерности numpy массивов и в объёмах вычислений.


    Для демонстрации простого, но максимально близкого к геологическому моделированию примера, на рисунке 6 представлена простая «пятискважинная» псевдо-модель.



    Рисунок 6. «Пятискважинная» псевдо-модель


    В определенных ячейках рассматриваемого объёма расположены «скважинные данные». Цветом показано распределение некоторого числового параметра, например пористости или проницаемости. Требуется провести моделирование неизвестных ячеек («межскважинного пространства») методами кригинга (рис. 7) и стохастического гауссова моделирования (рис. 8).



    Рисунок 7. Применение метода кригинга на «пятискважинной» псевдо-модели



    Рисунок 8. Применение метода стохастического гауссова моделирования на «пятискважинной» псевдо-модели


    Как видно по рисункам 7 и 8 оба метода дают «похожую» картину моделирования, однако моделирование кригингом даёт чёткий детерминированный образ, который будет воспроизводиться каждый раз при последующем моделировании, а стохастическая реализация более «реалистична» и каждый раз уникальна.


    Пример из «реальной» практики


    Но самое интересное происходит именно на данных, приближенных к реальным объектам. Например, модель «X» абстрактного месторождения «Y», на которой цветом показано уже распределение некоторого геологического параметра в объёме конкретного пласта. Моделирование методом кригинга представлено на рисунке 9.



    Рисунок 9. Применение метода кригинга на модели "X"


    Области с синими и красными цветами связаны со скважинами, имеющими минимальные и максимальные значения параметра. Кригинг соединил скважины друг с другом плавными переходами без резких скачков, а пространство в областях, далёких от скважин — устремил к среднему значению (зелёный цвет). Но картина со стохастическим моделированием имеет немного другой характер (см. рис. 10).



    Рисунок 10. Применение метода стохастического гауссова моделирования на модели "X"


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


    Выводы


    Целью нашей статьи было показать, что в части пространственной интерполяции в геологии нет ничего принципиально сложного. Прототипы этих геологических операций часто умещаются на одном экране. И между тем, это действительно математическое сердце дорогих и сложных программных продуктов. Современное геологическое моделирование не ушло далеко от простейших линейных моделей, и здесь ещё много можно сделать.


    Разумеется, это не все интерполяционные модели, используемые в нефтяной геологии, мы не обсудили вопросы моделирования дискретных (категориальных) полей, негауссовых случайных полей и т. д. За кадром так же остались вопросы построения структурных моделей (геометрического моделирования нефтяных пластов), гидродинамическое моделирование и многое другое. Так что соответствующие статьи ещё появятся в нашем блоге. Хотелось бы получить от читателя отзыв, не был ли математический материал слишком сложным или непонятным? Совсем отказаться от математики не кажется правильным, потому что, по крайней мере, в нефтяном инжиниринге вырваться в конкуренции с существующими компаниями вперёд можно только за счёт синтеза следующих направлений:


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

    Поэтому последующие статьи на тему, как что-то узнать, не вставая с дивана, будут содержать экскурсы в решения математических и физических задач, и, возможно, вызовут у читателей дополнительный интерес к фундаментальным областям науки, без которых решение многочисленных задач разведки и добычи не представляется возможным.

    Комментарии 45

      +1
      Было сложно, но я подучил!
        0

        Спасибо, упрощено до предела и даже дальше, но смысл есть. Имхо, еще стоило бы отметить, что о точных значениях речь все же не идет:


        интерполяция в точку с уже известным значением должна давать это самое значение.

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


        Также можно продемонстрировать, как разительно отличаются модели разных масштабов (уровней генерализации):
        image
        Как ни удивительно, это одна и та же модель объемной плотности на разных масштабах — и да, в геологии исходные данные разномасштабные и это необходимо учитывать. Зачастую, данные используются «как есть» без учета масштаба, просто изменяя разрешение данных, с соответствующим результатом… Кригинг — один из методов осмысленной передискретизации разномасштабных данных, доступный геологам, хотя еще лучше будет выполнять обработку в спектральной области (хотя и здесь есть свои сложности — скажем, оконное Фурье преобразование имеет право на жизнь, но не классическое Фурье преобразование, и так далее).

          +3
          Сильно напрягался, но так и не понял, как ваши картинки относятся к теме статьи?
            +2
            Одна из основных задач в этой области — распространить данные на объем, зная только то, что содержится в скважинах. Это можно сделать очень большим количеством методов, и в частности, надо учитывать «степень» изменчивости свойства. Крайгинг пытается это делать математически и не всегда корректно с геологической точки зрения. Можно помочь алгоритму с изменчивостью другими способами, например задав поверхность- или свойство-тренд.
              +2
              Это все правильно и понятно, вопрос был, где в картинках из комментария уважаемого N-Cube «hard data» (фактические данные)? И если ее тут нет, то какое отношение данные картинки имеют к задаче интерполяции и статье о распространении свойств? image
              –1

              Как я написал выше, это иллюстрация к тому, зачем вообще придуман и нужен кригинг — разница в данных разных масштабов не просто есть, она огромна. Если вы посчитаете среднюю плотность одного и того же объема на кубиках с картинки, она очень отличается. Кроме того, вычисление среднего при размере стороны ячейки меньше 400м на левом кубике и меньше 6000м на правом не имеет физического смысла. В статье этого нет, как нет и того, каков должен быть размер ядра функции (просто безразмерная константа «20» указана в формуле ковариации, для пиксельных вычисления).

                +1
                Как я написал выше, это иллюстрация к тому, зачем вообще придуман и нужен кригинг — разница в данных разных масштабов не просто есть, она огромна.
                Что-то совсем непонятно, кригинг просто один из методов интерполяции, при чем тут разница в данных разных масштабов. Где в ваших картинках фактические данные и где интерполированные?

                Кроме того, вычисление среднего при размере стороны ячейки меньше 400м на левом кубике и меньше 6000м на правом не имеет физического смысла. В статье этого нет, как нет и того, каков должен быть размер ядра функции (просто безразмерная константа «20» указана в формуле ковариации, для пиксельных вычисления)
                Ну в стате и много другого нет, так и что? Кажется, автор не претендует на исчерпывающее руководство по написанию пакета геостатистического моделирования. Какую кстати размерность вы хотите для коэффициента 20 в ковариации? Это же расстояние в палеопространстве, а так как в геологии оператором перехода из современного пространства в палео является сетка, то это расстояние в количестве ячеек.
                  –1

                  На картинке модель геологической плотности в двух разных масштабах, вычисленная из одних и тех же исходных данных методом инверсии (пространственный спектр рельефа — дополнение спектра грубых гравиметрических данных — вычисление источников гравитационного потенциала). В том и дело, что исходные данные одни и те же и разный результат. Ну, как с длиной береговой линии.


                  Это же расстояние в палеопространстве, а так как в геологии оператором перехода из современного пространства в палео является сетка...

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

                    +1
                    На картинке модель геологической плотности в двух разных масштабах
                    Не понимаю, не бывает «геологической плотности», бывает плотность вещества или энергии или еще чего-либо, что у вас? Что значит в разных масштабах, не одинаковы физические размеры, тогда как вы их, вообще, сравниваете? Честно говоря, больше походит на случайный набор слов.

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

                    1) Построить геологическую сетку, которая фактически является оператором перехода из настоящего пространства (где всё меряется метрами, футами и др.) с разломами, эррозией, выклиниванием, сжатиями и растяжениями в палео пространство, т.е. во время когда эти породы образовывались. В палео пространстве ещё не было никаких разломов и других факторов изменения пластов со временем.

                    2) В пиксельной (как вы её называете) системе палеопространства, распространить (интерполировать) свойства известные в некоторых точках (hard data) на всю заданную модель.

                    Переход в палеопространство необходим, так как вся теория геостатистики строится на знании процесса осадконакопления.
                      –1

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

                        +1
                        Дык вы ответите на это?
                        Не понимаю, не бывает «геологической плотности», бывает плотность вещества или энергии или еще чего-либо, что у вас? Что значит в разных масштабах, не одинаковы физические размеры, тогда как вы их, вообще, сравниваете?


                        Еще раз, алгоритмы геологической интерполяции (базируются на предположении, что мы привели наше текущее пространство в его исходное состояние назад во времени (палеопространство), именно тогда еще не было разломов, других нарушений. Именно в палеопространстве определяются радиусы вариограмм и др. параметры.
                        В палеопространстве можно использовать размерные величины, но только смысла в этом нет (так будет сложнее и дольше, но не будет точнее), так как обезразмеривание можно сделать на этапе перехода от современного в палеопространство, для чего используется структурная геологическая модель и геологическая сетка.
                          –1

                          Геологическая плотность — плотность при заданном масштабе геологического исследования. В статье же сказано про зависимость от масштаба.
                          Про палеопространство вы полную ерунду пишете, прикрывая невежество минусами ответам на свои посты. Нет никакого безразмерного идеального палеопространства без нарушений — нарушения были и есть всегда, а палеопространство это просто состояние в какой-то давний момент времени. Вы в самом деле считаете, что миллион-десять-сто назад длина была безразмерной? Безразмерная только вычислительная модель, а не физический мир.

            +1
            Конечно, выше подразумевается доверительный интервал

            Насколько я понял, подразумевается именно то, что написано: интерполятор в точках известных значений должен выдавать именно известные значения, без какого-либо доверительного интервала.
              –2

              Есть ошибки (малого количества) измерений керна и (множества) измерений геофизических параметров скважины по глубине (каротаж), по которым потом с определенной погрешностью вычисляются значения параметров для модели, и ко всему этому добавляется ошибка модели. Ни о каких точных значениях речь не идет — и при моделировании необходимо учитывать ошибки исходных данных и оценивать ошибку моделирования. К примеру, ошибка в 10% в исходных данных может привести к кратной разнице в результатах и т.п. Более того, полная ошибка в результатах модели может оказаться больше полученного результата, что делает его и вовсе бессмысленным. Это большая разница между инженерными дисциплинами, оперирующими данными реального мира, и программированием, оперирующим идеальными (ну почти, с учетом представления числа) значениями.

                +1
                Я, честно говоря, перестал вас понимать. Есть утверждение о том, что процедура интерполяции всего поля в точках известных нам должна воспроизводить известные нам значения. Вы не согласны с тем, что кригинг или метод обратных расстояний это делает? Или вы не согласны с таким требованием вообще и считаете его излишним?
                  –2

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

                    +1
                    Давайте я попробую ещё раз, последний. Вы берете цитату из статьи:
                    интерполяция в точку с уже известным значением должна давать это самое значение

                    Отвечаете на неё:
                    конечно, выше подразумевается доверительный интервал

                    Рискну ответить за автора (вы ведь тоже беретесь утверждать, что подразумевает автор!): никакого доверительного интервала не подразумевается. Упомянутые в статье интерполяторы обязаны воспроизводить значения hard-data точно, никак не интерпретируя погрешности приборов и так далее. Можно узнать ваш ответ на два простых вопроса, которые я упомянул сообщением выше?

                    Поэтому и множество стохастических моделей

                    Представьте себе гипотетический вариант, при котором используемые приборы имеют погрешность, которой можно принебречь. Раз вы пишете «поэтому» — получается, тогда множество стохастических моделей и не нужно использовать?
                    Конечно нет! Множество стохастических моделей нужно использовать не потому, что приборы врут, а потому что кригинг сглаживает!
                      –1
                      Упомянутые в статье интерполяторы обязаны воспроизводить значения hard-data точно

                      Нет, конечно, это не так. Объяснять нужно или сами разберетесь?


                      никак не интерпретируя погрешности приборов и так далее.

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


                      не потому, что приборы врут, а потому что кригинг сглаживает!

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

                        0
                        Кригинг и метод обратных квадратов расстояний являются точными интерполяторами? Они воспроизводят значения в тех точках, по которым построены?
                          –1

                          Вы отвечаете вопросом на вопрос. Но мне не сложно пояснить — метод обратных квадратов (как и других степеней, вообще говоря, посмотрите, что такое IDW) не определен в точке, где задано исходное значение. Знаете, на ноль делить нельзя… что-то забрезжило в сознании? Так вот, обычно из нерегулярной сетки пересчитывают в регулярную, так что значения исходные и интерполируемые пространственно не совпадают. Также, этот метод можно использовать для сглаживания на регулярной сетки — выкалывая центральную точку, то есть совпадающая исходная точка игнорируется. В результате получается сглаженная поверхность с ограниченным пространственным спектром. Итого, исходные и результирующие данные чаще всего определены в разных точках, но даже если они определены в одних и тех же точках, они не совпадают. Попытка взять совпадающие исходные точки «как есть» приведет к нарушениям гладкости поверхности и бесконечному спектру. Посмотрите спектры триангуляционной поверхности и интерполированной методом обратных квадратов, к примеру. Теперь понятно?

                            +1
                            Алексей! Раз уж мы на более личное общение переходим… Вы спрашиваете «теперь понятно?», но по факту делаете всё, чтобы было не понятно. У вас и у моих коллег, которые эту статью написали — противоположные цели. Вы, видимо, фрилансер, и поэтому намеренно всё усложняете, набивая себе цену; показываете, как всё сложно, и как с этой сложностью умеете справляться именно вы. Мои коллеги написали эту статью наоборот, чтобы показать, что в основе всего этого лежат простые вещи. Я как преподаватель, полностью разделяю этот подход — сначала надо разобраться в простых вещах, а только потом переходить на более сложные.

                            Есть простая математическая постановка задачи интерполяции. Есть конкретный математический метод: метод обратных квадратов расстояний. Да, он принудительно в точке имеющихся значений присваивает заданное значение, и об этом в статье написано — но ведь и предел при приближении к точке сходится туда же! Есть конкретный математический метод: кригинг. Он потому и называется точным интерполятором, что воспроизводит переданные в него значения точно. А вы вместо этого начинаете разводить наукообразную демагогию на тему того, как всё на самом деле сложно. Мы знаем, что всё сложно — но пытаемся объяснить простые вещи, которые лежат в основе!
                              –2

                              Вы бы воздержались от оскорблений, тем более, раз вы представитель компании. Тем более, что я статью не критиковал — по себе знаю, как сложно объяснить все это простым языком. Если уж в мой профиль заглянули, то в геологии я разрабатываю и публикую открытое ПО и модели, см. мой гитхаб: https://github.com/mobigroup Так что я на вашу должность не претендую, если что:) А если профиль весь прочитаете, то и вовсе станет стыдно хвалиться, что вы где-то преподаете.
                              Далее, далеко не все интерполяции в заданных точках дают исходные значения. Выше я написал про метод обратных расстояний с выколотой средней точкой — который по определению дает отличающееся значение в исходных точках. Зачем это нужно, я тоже написал. Более того, приведенные в статье интерполяции из неравномерно заданных значений на равномерную сетку также дают сглаженные поверхности — что и сказано в статье, при этом, оставлено за кадром, какие проблемы возникают, если одна или несколько исходных нерегулярных точек совпадут с точками регулярной сетки. При этом, ценны методы получения исходных значений в исходно заданных точках и при этом гладкой поверхности, и вот здесь речь заходит про криггинг. И об этом тоже сказано в статье. С чем вы спорите?


                              P.S. Почему я вообще трачу на эту дискуссию время, так это потому, что периодически есть необходимость просто и понятно описывать разные алгоритмы. Но конкретно в блоге данной компании — токсичность удивляет.

                                +1
                                Не я первый перешел на язык типа «что-то забрезжило в сознании», я пытаюсь подстраиваться под стиль того, с кем общаюсь. Но с вами общаться сложно. Когда вам задаёшь конкретный вопрос, который вроде бы должен читателям что-то прояснить, вы уходите в сторону и пускаетесь в усложняющие рассуждения.

                                Где-то в статье сказано, что все интерполяторы в заданных точках дают исходные значения? В статье сказано, что естественно это потребовать, с прицелом на те интерполяторы, которые рассматриваются ниже. Метод обратных расстояний имеет выколотую точку, но чему равен предел при стремлении к выколотой точке? А кригинг, я всё-таки хочу у вас узнать, даёт-таки исходные значения в исходных точках?

                                P.S. К сожалению, если посмотреть со стороны, то токсичность появляется у статей в этом блоге только в комментариях, и по удивительному совпадению, в ваших.

                                  –1

                                  Не вижу никаких усложняющих рассуждений. Давайте на примере. Даны два значения плотности — 2.95 и 3.05, точность определения каждого плюс-минус 10%. Так значение 3.0 для этих двух величин — это интерполяция или аппроксимация? По вашему мнению, это значение совпадает или не совпадает с измеренными данными? Сколько можно сделать интерполяций, при наличии разброса значений?


                                  Если непонятно, откуда ошибка, поясню еще раз. Пусть 5% ошибка определения плотности по несколько поврежденным образцам керна (см. статью), еще 5% ошибки добавляется из-за вычисления геологических параметров типа плотности по геофизическим измерениям (каротажным диаграммм) проводимости и проч. для диапазонов и скважин, где анализов керна нет (см. статью). Это еще скромная ошибка, на самом деле.


                                  В статье речь про интерполяцию кригингом, так что вычисляются только значения между измеренными. Но, как я уже писал выше, это не важно — при переходе от нерегулярной сетки к регулярной все значения на регулярной сетке являются интерполированными, то есть исходные значения не попадают в итоговые. Дело не в этом. Кригинг это такая интерполяция, которая минимизирует ошибку интерполированных данных в предположении гауусовости пространственного распределения значений. Таким образом, и исходные значения, заданные доверительным интервалом, прекрасно подходят, а «магический» коэффициент в ядре криггинга из примера статьи (20) и есть мера разброса значений! Даже в вики об этом явно сказано:
                                  image


                                  The kriging interpolation, shown in red, runs along the means of the normally distributed confidence intervals shown in gray.

                                  Но что делать, если распределение значений не гауссово? А если оно гауссово, но зависит не от расстояния? А если оно должно быть гауссовым, но при этом есть разрывы непрерывности? Вот тут и пригодится стохастическое моделирование, где можно и различные виды законов распределения попробовать, и разрывы непрерывности задать и так далее. И, разумеется, метод Монте-Карло намного более знаком каждому, нежели термин стохастического моделирования. И вы еще говорите, что это я что-то усложняю:)

                                    0
                                    Так значение 3.0 для этих двух величин 2.95 и 3.05 — это интерполяция или аппроксимация?

                                    По моему мнению это ни то ни другое, это просто третье число. По определению, оно не совпадает ни с тем, ни с другим, потому что число X по определению совпадает только с самим собой. Значение — это и есть значение, то есть просто число, оно совпадает только с самим собой. Интерполяций на двух ваших числах 2.95 и 3.05 нельзя сделать ни одной, потому что для интерполяции функции нужна независимая переменная и зависимая, а вы мне сообщили почему-то только видимо два значения зависимой переменной.

                                    Теперь мои вопросы:
                                    1) Какое отношение кригинг имеет к регулярной и нерегулярной сетке?
                                    2) Почему на вашей картинке серые доверительные интервалы схлопываются в точку там, где по оси X расположены красные квадраты, и почему схлопнутая точка доверительного интервала нулевой длины магическим образом совпадает с красным квадратом?
                                      –1

                                      Это значения с заданной погрешностью, то есть интервалы. Вы же любите безразмерные сетки — ну и определите сетку так, что эти значения находятся в соседних ее узлах.
                                      1) интерполяция нужна для получения приблизительных значений там, где они не определены исходно — то есть в узлах сетки для вычислений модели, которая, как правило, выбирается регулярной. Скажем, нужно вам моделировать миграцию нефти — необходимо пористость и прочие параметры в каждом узле расчетной сетки задать.
                                      2) картинка из википедии, статья кригинг, об этом явно сказано выше. Вам непонятно, что в вики написано или вы не смотрели?

                                        +1
                                        1) Вы меня с кем-то путаете, и автора статьи тоже. Где в формулах в статье хоть что-то сказано про сетку? Разве где-то в статье говорится о переходе от точечных замеров к сетке? Что для ваших двух значений является независимой переменной?
                                        2) Почему на вашей картинке серые доверительные интервалы схлопываются в точку там, где по оси X расположены красные квадраты, и почему схлопнутая точка доверительного интервала нулевой длины магическим образом совпадает с красным квадратом?
                              +1
                              Но мне не сложно пояснить — метод обратных квадратов (как и других степеней, вообще говоря, посмотрите, что такое IDW) не определен в точке, где задано исходное значение.
                              Что-то не очень у вас с математикой, с чего вы взяли, что значение в исходной точке не определено?

                              IMG-9851
                                –1

                                Значение 1/R не определено при R=0. Поэтому в первоначальном определении кригинга значения в исходных точках доопределяются отдельно.


                                У вас же на картинке α (слева) приравнено к αk (справа), что не имеет смысла. На самом же деле, это просто тождество и слева тоже должно стоять αk, но какой в этом смысл? Домножая числитель и знаменатель на квадрат Rk, вы просто заменили 1/Rk на коэффициент αk, пропорциональный 1/Rk, который вы не можете вычислить.


                                А кроме того, если считать по одной формуле и значения в исходных точках, доопределив соответствующие веса как единичные (вычислить их не удастся), из-за не нулевого вклада от остальных точек вы получите отличающиеся значения в исходных точках — против чего вы очень протестовали выше:)

                                  0
                                  Вообще взяв первую формулу, можно и числитель и знаменатель умножить на произведение всех Ri в квадрате, после чего вы получите эквивалентную формулу, у которой нет особых точек.
                                    +1
                                    Так, мне, наверное, надо объяснить, что альфа_i — это значения в известных точках i. альфа функция интерполяции методом обратных квадратов расстояний. Соответственно R_i — расстояние от искомой точки до известной точки i. Пока это просто формула. Умножим и числитель и знаменатель на R_k в квадрате, выражение 3, от этого же формула не поменяется. Внесем R_k в квадрате в сумму получим выражение 4. Теперь если будем рассматривать значение нашего интерполятора в точке k, то R_k будет равен нулю, а значит в выражении 4 суммы и в числителе и в знаменателе будут равны 0, а значение в точке k получим равным альфа_k, что и требовалось доказать.
                                    З.Ы. при чем тут кригинг?
                                      –2

                                      Вы зачем-то пытаетесь доказать тождество альфа_k равно альфа_k. Так чему равно это самое альфа_k? А равно оно 1/(R_k^2) и вычислить его вы не можете. И вы еще заявляете, что это у кого-то другого плохо с математикой?:)

                                        0
                                        Не совсем понимаю предмет дискуссии. В методе обратных расстояний (формула 2 в статье) в точке k очевидно значение не определено, но предел посчитать легко и он равен zk. Если вы собираетесь с этим спорить, мы не будем поддерживать больше эту дискуссию. И это не кригинг, вы наверное выше описались.

                                        Возникло несколько слов про доверительные интервалы. В этих математических моделях нет никаких доверительных интервалов и формулы точны. Если рассматривается определённая погрешность стоящая за измерениями, которые мы здесь интерполируем — эти вопросы нужно решать отдельно и не смешивать. Тоже самое относится к разломам.
                                          –1

                                          То есть выдаваемое выше за доказательство тождество от olegbor вы подменяете на формулу 2 из статьи? Просто отлично. Посмотрите выше, что в этой ветке обсуждается.


                                          Да, верно, про кригинг это другая ветка дискуссии, каюсь, попутал.

                                  0
                                  метод обратных квадратов (как и других степеней, вообще говоря, посмотрите, что такое IDW) не определен в точке, где задано исходное значение


                                  Давайте я сделаю вид, что ничего в этой теме не понимаю, и схожу хотя бы в вики на тему IDW. Что я там увижу?

                                  image

                                  Как же так, вы же говорите, что он не определён в точке? Скажете «не определён, а доопределён»?
                                    –2

                                    Выше в этой же ветке я так и сказал:


                                    доопределив соответствующие веса как единичные

                                    В ответ на olegbor:


                                    Что-то не очень у вас с математикой, с чего вы взяли, что значение в исходной точке не определено?

                                    Если что, «доопределить» это математический термин, вовсе не идентичный термину «определить».
                                    P.S. Вы там все такие, любой ценой стараетесь запутать и оскорбить собеседников?

                                      0
                                      Алексей! Я привел формулу из википедии, куда уж примитивнее. В этой формуле интерполяционная функция методом IDW определяется именно так, двояко. Еще раз, по вики, IDW определяется как ОБЕ эти формулы вместе с where. Вы говорите, как я и подозревал, что первое условие — это «определение», а второе условие — это «доопределение», по-вашему «вовсе не идентичное». Интересно, если бы там три было условия, как бы вы третье назвали, «додоопределение»?

                                      P.S. Мы все эти статьи пишем, чтобы попытаться систематизировать, упростить и объяснить. Но во все эти статьи вы приходите с комментариями, суть которых в том, чтобы показать, что всё гораздо сложнее — то есть действуете в противоположном нашему направлении. Я уже говорил, понятно, что вам, как фрилансеру, нужно себя зарекомендовать. Может быть уже хватит? Оставьте нас уже с нашими заблуждениями и примитивными представлениями? С удовольствием загляну в вашу статью на тему интерполяции.
                    –3

                    Кто бы то ни был (Я имею в виду автора статьи). Всё бы ничего, да вот одна проблема, текст выглядит вымученным, написанным через силу, но, якобы, с претензией на литературность. Нет. Нет её тут. Автор, видно, тебе привычно писать около научно-официозным стилем, так лучше бы им и написал, а не мучил себя, меня и читателей. Сама статья хорошая, тема интересная.

                      0
                      Спасибо за статью!
                      Было бы интересно еще услышать про объектное моделирование и нейронные сети — сейчас это модная тема и на многих геологических конференциях позиционируется как «наше все».
                        0

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

                          0
                          Автор рассказывает не о том, как данные для нейронных сетей готовить, а о том, как интерполировать по-разному можно неизвестные данные между известными точками.
                            –1

                            Не рассказывает, на самом деле. Потому что при наличии в реальном мире разрывов непрерывности все это не работает. И тот самый дорогой софт, о котором с таким пренебрежением автор то и дело упоминает, просто умеет работать с таким вещами, о которых автор, видимо, даже не догадывается.

                              0
                              Откуда взялся разрыв непрерывности? Разлом? До того, как разлом образовался, непрерывность была?
                                –1

                                Не было, разумеется, поскольку никогда литосфера не была статичной. Древние разломы перекрываются современными, геологические блоки двигаются, литосферные плиты раздвигаются или сдвигаются и трутся и т.п. Постройте интерферограммы любого участка земной поверхности за полвека (за это время можно найти данные) и вы увидите, что изменения происходят непрерывно (в том числе, очень значительные изменения — сдвиги на десятки и сотни метров вполне возможны).
                                Кстати, немало публикаций на тему выделения палеоразломов — это важная и нужная задача, а палеоповерхности выделить проще. Вот пример из практики — в Африке не редкость газ, прямо выходящий на поверхность, при этом горизонты залегания привязаны к палеоповерхностям, а точки выхода на поверхность к современным разломам. Притом, этот газ мигрировал из кристаллического фундамента по палеоразломам, уже «заросшим», в палеоловушки, потом они нарушились более современными разломами...

                                  –2

                                  Вот, как раз в тему статья про кригинг интерполяцию при наличии разломов, попалась по тегу VTK, как ни странно:
                                  3D Visualization of Stratum with Faults Based on VTK
                                  Тема-то общеизвестная и актуальная… в отличие от просто формулы из вики.

                        Только полноправные пользователи могут оставлять комментарии. Войдите, пожалуйста.

                        Самое читаемое