Мы уже обсуждали современные методы в геологии в статье Геология XXI века как наука данных о Земле на примере модели землетрясения в горном массиве Монте Кристо в Неваде, США 15 мая 2020 года магнитудой 6.5 баллов. И все бы хорошо в этой модели, да вот только самое интересное смещение геологических блоков и "дыхание гор" нам схематично указал опытный геолог. Самое же важное заключается в том, что современная вычислительная геология (включая геофизику, моделирование и визуализацию) позволяет создать динамическую (4D) геологическую модель и наяву увидеть происходящие геологические процессы.

Геологическая модель с интерферограммой на поверхности рельефа по данным радарной спутниковой съемки, где на шкале Density Anomaly,% является характеристикой неоднородности геологической плотности и черная сфера в центре указывает координаты эпицентра землетрясения, расположенного на глубине 2.8 км.

Поскольку в указанной выше статье мы уже рассмотрели статичную модель, сразу перейдем к динамической модели и ее визуализации. Как обычно, воспользуемся для этого Open Source программой ParaView и моим расширением для ГИС данных N-Cube ParaView plugin for 3D/4D GIS Data Visualization. Вот как выглядит проект ParaView:

Напомню, что геологическую модель мы создаем методом так называемой геофизической инверсии, когда на основе данных гравитационного поля на поверхности Земли вычисляем соответствующее распределение плотности под этой поверхностью. Увы, но измерения непосредственно гравитационного поля (или нужной нам вертикальной его компоненты) с такой точностью и периодичностью не производятся, поэтому воспользуемся заместо этого открыто доступными регулярными радарными снимками. Дело в том, что пространственные спектры гравитационного поля, рельефа и радарных (и оптических) снимков практически эквивалентны, что и дает возможность восстановить распределение плотности с точностью до множителя. Если вам интересны детали, то в GitHub репозитории GIS Snippets доступны Jupyter Python 3 ноутбуки с соответствующими моделями (и ссылками на теоретические основы). Спутниковая интерферограмма получена средствами замечательного открытого тулкита GMTSAR.

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

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

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

Красным цветом на интерферограмме показаны поднявшиеся участки (в данном случае, на 20-30 см), а синим опустившиеся (на 15-20 см). Смотрите подробнее в статье Общедоступные данные дистанционного зондирования Земли: как получить и использовать

В заключение скажу, что сам я очень рад увидеть такую динамическую модель несмотря на более чем десятилетнюю работу с геологами и десятки построенных 3D геологических моделей, мне все еще были не очень понятны их объяснения динамики смещения, поворотов и скольжения геологических блоков Еще замечу, что данная динамическая модель построена исключительно для удовлетворения моего научного интереса, посколько для геолога здесь и так все очевидно.

В заключение, приглашаю всех посетить GitHub репозитории с множеством геологических моделей и их визуализацией в Blender и ParaView, а также примерами анализа пространственных спектров, синтеза гравитационных данных высокого разрешения на основе данных дистанционного зондирования и другими вычислениями, в том числе, выполняемыми на геоиде средствами PostgreSQL/PostGIS. Также смотрите готовые визуализации на YouTube канале.

Сегодня вместо обсуждения геологических моделей мы посмотрим пример их программирования в среде Jupyter Notebook на языке Python 3 и с библиотеками Pandas, NumPy, SciPy, XArray, Dask Distributed, Numba, VTK, PyVista, Matplotlib. Это довольно простой ноутбук с поддержкой многопоточной работы и возможностью запуска локально и в кластере для обработки больших данных, отложенными вычислениями (ленивыми) и наглядной трехмерной визуализацией результатов. В самом деле, я постарался собрать разом целый набор сложных технических концепций и сделать их простыми. Для создания кластера на Amazon AWS смотрите скрипт AWS Init script for Jupyter Python GIS processing, предназначенный для единовременного создания набора инстансов и запуска планировщика ресурсов на главном инстансе.

Визуализация с помощью Visualization Toolkit(VTK) и PyVista это уже далеко не Matplotlib

Идея сделать такой пример возникла у меня давно, поскольку я регулярно занимаюсь разнообразными вычислительными задачами, в том числе для различных университетов и для геологоразведочной индустрии, и знаком очень близко с проблемами переносимости и поддерживаемости программ, а также проблемами работы с так называемыми большими данными (сотни гигабайт и терабайты) и визуализацией результатов. Так что само собой появилось желание сделать ноутбук-пример, в котором коротко и просто показать и красивую визуализацию и распараллеливание и ускорение кода Python и чтобы этот ноутбук можно было без изменений запустить как локально, так и на кластере. Все использованные библиотеки доступны уже много лет, но мало известны, или, как говорится, они остаются широко известными в узких кругах. Оставалось лишь найти подходящую задачку, на которой все это можно показать и это было, пожалуй, самым сложным ведь мне хотелось, чтобы пример получился достаточно осмысленным и полезным. И вот такая задача нашлась рассмотреть моделирование гравитационного поля на поверхности для заданной (синтетической в данном случае) модели плотности недр и некоторые последующие преобразования с вычислением фрактального индекса по компонентам пространственного спектра и кольцевого преобразования Радона, как его называют математики, или Хафа, согласно компьютерным наукам. Замечательно то, что с популярными библиотеками Python эти преобразования делаются буквально в несколько строчек кода, что особенно ценно для примера. Поскольку моделирование поля в каждой точке поверхности требует вычисления для всего трехмерного объема, мы будем обрабатывать гигантский объем данных. Для визуализации используем человеколюбивую обертку PyVista для библиотеки VTK Visualization Toolkit, потому что писать код для последней это путь истинных джедаев кто хочет лично в том убедиться, смотрите мой модуль к ParaView N-Cube ParaView plugin for 3D/4D GIS Data Visualization, написанный как раз на Python + VTK.

Теперь предлагаю проследовать по ссылке на страницу GitHub репозитория или сразу открыть ноутбук basic.ipynb Надеюсь, код достаточно просто читается, остановлюсь лишь на нескольких особенностях. Запускаемый в ноутбуке локальный кластер dask предназначен для работы на многоядерных компьютерах, а вот для работы в кластере потребуется настроить подключение к его планировщику. В упомянутом выше скрипте AWS Init script for Jupyter Python GIS processing есть соответствующие комментарии и ссылки. В коде мы используем векторизацию NumPy, то есть передаем сразу массивы, а не скаляры, при этом пользуемся тем, что XArray объекты предоставляют доступ к внутренним NumPy объектам (object.values). Код NumPy ускорить непросто, но с помощью Numba и для такого кода можно получить некоторый выигрыш в скорости исполнения (возможно, даже около 15%):

from numba import jit@jit(nopython=True, parallel=True)def delta_grav_vertical(delta_mass, x, y, z):    G=6.67408*1e-11    return -np.sum((100.*1000)*G*delta_mass*z/np.power(x**2 + y**2 + z**2, 1.5))

Для перебора всех точек на поверхности куба и вычисления для каждой такой точки гравитационного воздействия от каждой точки куба (это, кстати, получается пятикратный интеграл) мы пишем:

def forward_gravity(da):    (da_y, da_x, da_z) = xr.broadcast(da.y, da.x, da.z)    deltagrav = lambda x0, y0: delta_grav_vertical(da.values.ravel(), (da_x.values.ravel()-x0), (da_y.values.ravel()-y0), (da_z.values.ravel()-0))    gravity = xr.apply_ufunc(deltagrav, da_x.isel(z=0).chunk(50), da_y.isel(z=0).chunk(50), vectorize=True, dask='parallelized')    ...

Здесь xarray.broadcast с линеаризацией массивов функцией ravel() позволяют из трех одномерных координат x, y, z получить триплеты координат для каждой точки куба. Выражения da_x.isel(z=0) и da_y.isel(z=0) извлекают x, y координаты верхней поверхности куба, на которой и вычисляется гравитационное поле (точнее, его вертикальную компоненту, т.к. именно она измеряется при практических исследованиях и такие данные доступны для анализа). Функция xarray.apply_ufunc() весьма универсальная и одновременно обеспечивает векторизацию и поддержку параллельных ленивых вычислений dask для указанной коллбэк функции deltagrav. Хитрость заключается в том, что для выполнения вычислений на кубе для каждой точки поверхности нужно координаты поверхности передать в виде XArray массивов, а для использования dask они также должны быть dask массивами, что мы и обеспечиваем конструкциями da_x.isel(z=0).chunk(50) и da_y.isel(z=0).chunk(50), где 50 это размер блока по координатам x, y (подбирается в зависимости от размера массивов и количества доступных вычислительных потоков). Да, такая вот магия достаточно лишь использовать вызов chunk() для XArray массива, чтобы автоматически превратить его в dask массив.

Обратим внимание, что dask-вычисления по умолчанию являются ленивыми (отложенными), то есть вызов функции forward_gravity() завершается почти мгновенно, но возвращаемый результат является лишь оберткой, которая инициирует вычисления только при непосредственном обращении к данным или вызовом load(). При интерактивной работе это очень удобно, так как мы можем написать сложный пайплайн с большими наборами данных и для проверки и визуализации выбрать лишь маленький его кусочек, а при необходимости и запустить вычисления на полном наборе данных. К примеру, мне часто приходится работать с NetCDF датасетами глобального рельефа планеты и прочими в сотни гигабайт на своем ноутбуке визуализируя малую часть данных, а потом запускать уже готовый ноутбук в облаке для обработки всех данных. Таким образом, код для локальной работы и его продакшен версия ничем не отличаются. Главное, правильно настроить размеры dask блоков, иначе вся магия "сломается".

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

В заключение, приглашаю всех посетить GitHub репозитории с множеством геологических моделей и их визуализацией в Blender и ParaView, а также примерами различного анализа. Также смотрите готовые визуализации на YouTube канале.

Ранее в статьях мы уже обсудили доступные данные (результаты наземных и спутниковых гравитационных и магнитных измерений, ортофото и космические снимки, цифровые модели рельефа), теоретические подходы и методы обработки (интерферометрия, построение обратных геофизических моделей), обработку данных в ParaView (выделение изоповерхностей) и Blender (высококачественная визуализация и анимация подготовленных в ParaView данных) и даже посмотрели Python Jupyter notebook с вычислениями и визуализацией моделей (включая выделение изоповерхностей средствами библиотеки VTK). Осталось построенные геотермальные изоповерхности конвертировать в формат модели дополненной реальности и получить геотермальную модель в дополненной реальности(AR). Эта модель может быть легко просмотрена прямо из браузера на iOS/iPadOS или из загруженного файла на MacOS (поддержка AR добавлена два-три года назад, на старых устройствах для этого потребуется обновиться). Увы, стандартов много и удастся ли открыть модель на Windows или Android, я не знаю (напишите в комментариях, если можно добавить поддержку нужного стандарта каким-то софтом). Как всегда, исходная модель доступна на GitHub в репозитории ParaView-Blender в виде исходных STL/PLY файлов и проекта Blender.

AR Модель геотермального резервуара Лахендонг, полуостров Минахаса, Северный Сулавеси, Индонезия Замеры температуры по скважинам обозначены цветными дисками синим 0-150C (далеко от резервуара), белым 150-250C (переходная область вблизи от резервуара), красным 250-350C (внутри геотермального резервуара).

Вместо введения, или зачем все это нужно

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

К примеру, мы не можем пробурить скважину до центра Земли и собрать образцы, поэтому приходится придумывать другие способы исследования скажем, регистрировать волны мощных землетрясений и прочих геологических событий в разных местах планеты и сопоставлять, пытаясь определить, прошла ли волна напрямую (по хорде) через толщу земли или отклонилась, как волна замедлилась и так далее. В результате, когда была пробурена Кольская сверхглубокая скважина, вся наша устоявшаяся геологическая картина оказалась ошибочной неведомым образом температуры в глубине вдвое выше ожидаемых (к потоку тепла из глубинных слоев еще столько же добавляется в результате нагрева от радиоктивного распада), породы не плотные, а очень трещиноватые и заполненные водой (притом, на поверхности начали бурить сплошной и прочный Балтийский щит возрастом миллиарды лет), отсутствует ожидаемый переход между гранитами и базальтами на глубине около 7км (этот переход отлично виден на сейсмограммах и не было сомнений в его существовании), При этом, образцы пород (керн) удавалось получить в виде грязи или пыли, в зависимости от влагонасыщенности пород (то есть даже простое, на первый взгляд, измерение плотности пород задача весьма непростая). Почему я говорю именно про эту скважину? Потому, что это проект мирового масштаба и призванный (в свое время) показать всю силу советской науки, то есть выполнялся на "высшем уровне", в то время как результаты коммерческого бурения очень далеки от него по качеству сбора образцов и их анализа, и все подробности и неточности сбора и обработки данных не регистрируются.

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

Еще есть совпадения и ошибки интерпретации. Поскольку геологоразведка состоит из многих этапов (и никогда не заканчивается), зачастую интерпретацию данных каждого этапа используют на последующих без пересмотра и анализа всех собранных данных совместно. Рассмотрим на примере нефтяного участка в России, уже зарегистрированного с отличным ресурсным потенциалом (извлекаемыми запасами нефти). Дело было примерно так: разведка, как полагается, началась весной, пара пробуренных скважин оказались рабочими и показали отличное давление нефтяного пласта. Геологи проанализировали данные сейсмики и бурения и нарисовали модель мощной нефтяной линзы, посчитали извлекаемый объем (много!), с этими результатами зарегистрировали месторождение и участок был продан компании, занимающейся добычей нефти. Но вскоре что-то пошло не так и давление в скважинах упало. Не беда, главное, это отличное новое месторождение. У скважин попробовали менять режим работы, возможно, попробовали еще немного пробурить и так далее и на следующую весну они снова фонтанировали нефтью. Все бы хорошо, но через пару месяцев давление снова упало Проверка всех данных показала, что найденный "нефтяной пласт" не был проанализирован и оказался водным (что, кстати, подсказывал рельеф местности и водотоки), а нефть (более легкая, чем вода) находилась лишь в тонкой прослойке, и мало того, хорошее давление в пласте создавалось лишь весной за счет грунтовых вод (достаточно сравнить подъем поверхности по данным спутниковой интерферометрии, совпадающий с периодами высокого давления в скважинах). Почему не пробурили саму потенциальную нефтяную линзу для проверки? Все просто ее центр находится за пределами данного участка и принадлежит другому владельцу, так что там бурить не было возможности. Что ж, бывает, именно поэтому многие компании, занимающиеся разведкой полезных ископаемых, не занимаются их промышленной добычей.

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

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

Таким образом, из вышесказаного становится ясно, что никаким данным априори доверять нельзя, а нужно все проверять и сопоставлять, а при получении геологического результата снова перепроверять те данные, которые к нему привели. Именно поэтому нейросети и прочие недетерминированные методы применяются, в основном, для прикидки как правило, это далеко не лучшая идея засунуть в "черный ящик" мусорные данные ("мусор на входе мусор на выходе"). Спасают во всем этом калейдоскопе методы анализа и очистки данных (data science) и качественная визуализация, которой никогда не бывает много -на ней геолог видит намного больше, чем специалист по данным, который эту визуализацию и сделал. Вряд ли нам удастся выбрать идеальный ракурс для статичной модели с точки зрения геолога (притом, как вы уже догадались, три геолога попросят три разных ракурса). Возможности дополненой или виртуальной реальности позволяют легко и быстро масштабировать и перемещать модель, просматривать с любый ракурсов и делать скриншоты для их сверки с данными получается виртуальное геологическое исследование. Все это в общем-то понятно и очевидно, а действительно актуальным стало лишь с появлением широко распространенных стандартов, поддерживаемых во множестве смартфонов, которые сегодня есть почти у всех.

Технические подробности

Как говорится, все развивается по спирали, и с технологией виртуальной реальности с моделями формата VRML я столкнулся полтора десятка лет назад. Увы, это было практически бесполезно, поскольку ресурсов обычного десктопного компьютера не хватало на комфортную работу с моделью, а через несколько минут вынужденно медленного просмотра приходилось перезапускать программу просмотра, поскольку она выходила за пределы доступной оперативной памяти и дальнейшая работа становилась невозможной. Это уж не говоря о необходимости ставить какой-то сомнительный софт для просмотра таких моделей, причем программ просмотрщиков было много и они были несовместимы между собой Сегодня многое изменилось каждый более-менее современный MacOS компьютер и iOS/iPadOS устройство имеют поддержку просмотра моделей дополненной реальности, причем это очень удобно делать именно с мобильного устройства, которое при этом даже не грееется на ощупь и позволяет очень плавно и со всех сторон осматривать модель и с ней взаимодействовать (масштабировать, перемещать, вращать). Также возможно добавлять в модели различные триггеры событий и ссылки на веб-адреса, делать анимации и так далее. Все это и послужило причиной, почему я, после обновения на MacOS Catalina (мне и с предыдущей хорошо было, так что обычно я жду года после выхода, прежде чем обновляться на новую, уже стабильную систему) решил попробовать сделать такую модель в дополненной реальности. Чтобы упростить себе задачу, начал со статической модели, вдобавок, которая у меня уже готова в виде проекта Blender.

Apple предлагает набор средств разработки AR Creation Tools, из которых мне пока потребовался только консольный Python модуль USDZ Tools, а рекомендуемый Reality Composer потребовал установки среды разработки XCode (внимание: сразу после инсталляции занимает 30ГБ места на диске) и еще не пригодился. Отдельно устанавливаемый Reality Convertor умеет шуметь вентилятором и делает то же самое, что и Python модуль, а еще в нем можно красивый скриншот модели сделать (смотрите картинку в начале статьи).

Смотрите на GitHub модель и инструкции для переноса данных из ParaView в Blender и в AR модель (или сделанных с помощью библиотеки VTK в Jupyter notebook, как описано в моей предыдущей статье): ParaView-Blender В поддиректории "export" кратко приведены подробности переноса и подготовки данных. Сами исходные данные модели и проект Blender доступны в поддиректории "Minahasa". Там же доступны "сырые" данные в виде скриптов Google Earth Engine (GEE), GeoJSON, TIFF, NetCDF файлов обратной геофизической модели для исходного проекта ParaView. Вот как выглядит нужный нам резервуар:

Рендеринг из проекта Blender на GitHub для AR упростим проект

Формат экспорта gITF 2.0 позволяет разом перенести множество объектов из Blender, команда конвертации приведена на гитхабе по ссылке выше. При этом, необходимые текстуры требуется предварительно подготавливать и сохранять в отдельные файлы с помощью их "прожига" (Bake). Для переноса в формат модели AR пришлось обойтись диффузным шейдером и "прожигом" (Bake) только диффузного света согласно цветам узлов модели как показано ниже:

Проект Blender с готовой для экспорта в AR моделью за основу взят проект из репозитория ParaView-Blender

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

Заключение

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

В прошлой статье Google Earth Engine (GEE) как общедоступный суперкомпьютер речь шла про работу в облачном редакторе GEE, где для доступа достаточно лишь наличия Google почты. Если потребности ограничиваются разовыми задачами и гигабайтами извлекаемых данных, то этого вполне достаточно. Но для автоматизации множества даже мелких задач облачный редактор не лучший способ работы и, тем более, когда требуется многократно получать растры суммарным размером в терабайты. В таких случаях потребуются другие инструменты и сегодня мы рассмотрим возможности доступа из консольных shell и Python скриптов и Python Jupyter notebook.

На скриншоте Python Jupyter ноутбук, где растр с данными о плотности населения за 2020 год из каталога Earth Engine data Catalog: WorldPop Global Project Population Data отображен на карте OpenStreetMap

Введение

Как говорится, аппетит приходит во время еды. Вот и я, хотя давно и с удовольствием работаю с Google Earth Engine (GEE), только недавно столкнулся с задачами, требующими регулярного извлечения терабайт данных для их последующей обработки. К счастью, это вполне реально. Более того, на GEE все нужные данные представлены в одном месте и регулярно обновляются, что делает его оптимальным источником. Конечно, мне понятно, насколько мало специалистов работают с такими наборами пространственной информации, и уж точно они не ищут статьи на русском языке (потому, что их нет). С другой стороны, специалисты по машинному обучению (ML) часто жалуются на недостаток данных, так вот же вам настоящий Клондайк! Есть разные варианты, как реализовать машинное обучение на этих данных можно средствами GEE, можно обратиться к функциям Compute Engine или делать все самому. Впрочем, это отдельная большая тема, так что мы ограничимся лишь получением данных.

Облачные акаунты Google

Для работы нам понадобится подключение к облачным аккаунтам Google Cloud SDK и установленный консольный google-cloud-sdk. При наличии гугл почты у нас уже есть один (или несколько) персональных аккаунтов. Просмотреть список доступных аккаунтов и переключаться между ними можно с помощью консольной команды:

$ gcloud auth listCredentialed accounts: - youremail@gmail.com (active)To set the active account, run $ gcloud config set account <account>

Здесь уже отображена подсказка для переключения между аккаунтами:

$ gcloud config set account <account>

После переключения нужно авторизоваться в аккаунте :

$ gcloud auth login

Доступ к облачным хранилищам

Обычных почтовых аккаунтов достаточно для доступа к облачным хранилищам buckets и Google Drive, куда можно сохранять данные GEE из облачного редактора GEE. Заметим, что при доступе через API данные возможно сохранить сразу локально.

Cохранение данных GEE на buckets выполняется с помощью функций Export.table.toCloudStorage и Export.image.toCloudStorage и используется в случаях, когда планируется дальнейшая работа с файлами в облаке Google Compute Engine. Управлять этими файлами можно с помощью утилиты gsutil, например:

$ gsutil du -h gs://gcp-pdp-osm-dev-earth-engine

Эта команда покажет список файлов и их размеры в удобных человеку единицах (см. ключ -h). С помощью утилиты gsutil можно осуществлять различные файловые операции, ключи достаточно интуитивны (cp, rm,...), подробности смотрите в справке указанной утилиты.

Для сохранения из GEE на Google Drive доступны команды Export.table.toDrive и Export.image.toDrive, дальнейшее управление файлами доступно в веб-интерфейсе или с помощью дополнительных утилит и приложений. Преимущественно сохранение на Google Drive используется для скачивания готовых файлов.

Доступ к GEE через API

Для не интерактивной работы с Google Earth Engine (GEE) рекомендуется создать так называемый сервисный аккаунт вида my-service-account@...iam.gserviceaccount.com: Create and register a service account to use Earth Engine. Для доступа к GEE требуется в свойствах созданного аккаунта зайти на вкладку KEYS и создать и скачать JSON файл его ключа, а также зарегистрировать аккаунт на странице Register a new service account. Теперь с этим ключом можно авторизоваться с помощью следующего Python кода:

import eeservice_account = 'my-service-account@...iam.gserviceaccount.com'credentials = ee.ServiceAccountCredentials(service_account, 'privatekey.json')ee.Initialize(credentials)

вместо использования функции Python API ee.Authenticate() с интерактивной авторизацией каждого сеанса. Смотрите также консольную авторизацию

$ earthengine earthengine --ee_config

Аналогично, ключ сервисного аккаунта позволяет авторизоваться и при использовании Python GDAL:

import osfrom osgeo import gdalos.environ["GOOGLE_APPLICATION_CREDENTIALS"] = "my-service-account.json"

Или консольных утилит GDAL:

export GOOGLE_APPLICATION_CREDENTIALS=my-service-account.json

Извлечение растров GEE

Для получения растров мы будем говорить о методе API Method: projects.assets.getPixels Такой доступ обеспечивает высокую скорость передачи данных и позволяет копировать огромные объемы данных, но отдельными блоками размером не более 32MB. К счастью, проект GDAL уже предоставляет обертку для этого API, так что достаточно одного вызова нужной утилиты или функции для передачи всего растра целиком.

Возьмем практический пример для консольных утилит GDAL и на Python. Просмотрим набор данных WorldPop/GP/100m/pop и извлечем один из найденных растров за 2020 год. Растры в этом наборе варьируются в размере по порядку величины от мегабайт до гигабайт, для примера выберем один из небольших:

export GOOGLE_APPLICATION_CREDENTIALS=my-service-account.json# fetch collectionogrinfo -ro -al "EEDA:" -oo COLLECTION=projects/earthengine-public/assets/WorldPop/GP/100m/pop -where "year=2020" # show one raster infogdalinfo "EEDAI:projects/earthengine-public/assets/WorldPop/GP/100m/pop/ZWE_2020"# fetch one raster to local drivegdal_translate "EEDAI:projects/earthengine-public/assets/WorldPop/GP/100m/pop/ZWE_2020" ZWE_2020.tif

Аналогично на Python:

import osfrom osgeo import ogr, gdal# define service account keyos.environ["GOOGLE_APPLICATION_CREDENTIALS"] = "my-service-account.json"# fetch collectiondriver = ogr.GetDriverByName('EEDA')ds = driver.Open('EEDA:projects/earthengine-public/assets/WorldPop/GP/100m/pop')layer = ds.GetLayer()# filter collection by attributelayer.SetAttributeFilter('year=2020')# select 1st rasterfor feature in layer:    name = feature.GetField("name")    crs = feature.GetField("band_crs")    print ('raster name and crs:',name, crs)    break# fetch 1st raster from the collection to arrayds = gdal.Open(f'EEDAI:{name}')band = ds.GetRasterBand(1)array = band.ReadAsArray()print ('raster shape:', array.shape)

Заключение

Выше мы рассмотрели достаточно продвинутые техники работы с Google Earth Engine. Если вы обращаетесь к GEE регулярно, намного удобнее использовать шелл скрипты и Python Jupyter ноутбуки и иметь возможность сохранять практически произвольные объемы данных локально без промежуточных облачных хранилищ и без ожидания выполнения очереди заданий экспорта, которое может затянуться. Отдельно отмечу, что вовсе не обязательно извлекать сырые данные можно их предварительно обработать средствами GEE. За подробностями серверной обработки, работы с GDAL и отображения данных обратитесь к ссылкам ниже.

Мне было бы интересно получить обратную связь от читателей: стоит ли обращаться к более сложным темам или и это уже за гранью того, что интересует русскоязычную аудиторию? Знаю, что немало читателей здесь используют Google Transtale и подобные переводчики, возможно, стоит сразу писать на английском на LinkedIn, как я уже делаю с публикациями по геофизике.

Ссылки

EEDAI Google Earth Engine Data API Image

Raster API tutorial

Vector Layers

Using GDAL / OGR for Data Processing and Analysis

Raster Layers

Raster (gridded) dataset handling

FROM GEE TO NUMPY TO GEOTIFF

How to load GeoJSON files into BigQuery GIS

В предыдущих статьях Google Earth Engine (GEE) как общедоступный суперкомпьютер и Google Earth Engine (GEE) как общедоступный каталог больших геоданных мы познакомились со способами удобного и быстрого доступа к каталогу космических снимков и их обработки. Теперь мы можем искать питьевую воду, различные минералы и вообще много всего. А еще можем вооружиться методами машинного обучения (ML) и сделать свою собственную карту сокровищ прогноз для поиска золотых месторождений в любом месте мира. Как всегда, смотрите код и исходные данные (синтетические, конечно, ведь реальные данные буквально на вес золота!) на GitHub: AU Prediction (ML)

На острове Западная Сумбава с помощью построенного классификатора выделены прогнозируемые золотоносные участки.

Введение

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

Вручную (хоть и на компьютере, оксюморон) все это делается, но делается медленно и разными несогласованными методами с несомненным преобладанием человеческого фактора. Хорошо нам знакомая и необходимая итеративная работа невозможна, так как нет автоматизированного способа повторять подготовку данных и их анализ. Для нас задача становится решаемой благодаря наличию доступа к единому многолетнему архиву всех необходимых космических снимков, рельефа и производных от него материалов, а также синхронизированных по времени масок облачности и классификаторов (вода, суша, растительность и других) на Google Earth Engine.

Подготовка данных

Обычно геологи работают с одним или несколькими (выбранными на свой вкус) снимками лишь одной спутниковой миссии, мы воспользуемся техническим преимуществом платформы Google Earth Engine (GEE) и из тысяч доступных космических снимков сделаем качественные композитные растры для разных миссий (Landsat 8, Sentinel 2, ...) и дополним их детальным цифровым рельефом местности и некоторыми его характеристиками.

Остановимся подробнее на том, зачем нужны снимки с разных миссий, хотя и разрешение и спектральные каналы у них более-менее одинаковые. На самом деле, нет разная оптика и ширина спектральных каналов и алгоритмы обработки обуславливают существенную разницу в результирующих изображениях. Ранее я уже публиковал свои результаты сравнения значений спектральных каналов разных миссий для классов поверхности (земля-вода-растительность) как раз для исследуемой территории: Compare Spectrograms of Hyperspectral and Multispectral Satellite Missions:

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

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

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

Построение классификатора и анализ территории

На гитхабе представлено два ноутбука, в первом из которых GEE_export.ipynb подготавливаются и скачиваются необходимые данные с Google Earth Engine, а во втором West Sumbawa AU 3-Class Prediction Synthetic.ipynb выполняется обучение классификатора на полученных растрах и их последующая классификация. В коде даны ссылки на все используемые датасеты GEE, так что, при желании, их несложно заменить на другие. Поскольку использованный линейный классификатор достаточно прост и требует подбора лишь одного параметра, этот подбор выполняется автоматически непосредственно в ноутбуке.

Отмечу, что желательно использовать физические характеристики поверхности функции от рельефа, отражательную способность для разных длин волн и так далее. Тем не менее, вместо отражающей способности земной поверхности (surface reflectance) для миссии Landsat-8 мы используем данные TOA (top of atmosphere), поскольку для этих спутников SR вычисляется по TOA недостаточно аккуратно и нужная нам корреляция теряется. Для Sentinel-2 такой проблемы нет.

Заключение

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

P.S. С днем Геолога!

Сегодня мы будем искать полезное ископаемое золото с помощью открыто доступных на платформе Google Earth Engine (GEE) данных, используя геологическое моделирование и последующую классификацию методом опорных векторов для предсказания золотоносных участков по построенной геологической модели. Нам понадобится рельеф ALOS разрешением 30 м, радарные снимки Sentinel-1 SAR разрешением 10 м и оптические снимки Sentinel-2 10 м (только для визуализации). Точность классификатора получилась равной 97.77% и, самое главное, результат соответствует ожиданиям геолога найденные участки на самом деле очень перспективны.

Красно-белым шариком отмечен участок для детального исследования

Введение

Геологически Индонезия находится в так называемом Тихоокеанском рудном поясе, освоение которого первыми, вероятно, начали инки, ацтеки и майя на американском континенте, а за ними эстафету приняли испанские конкистадоры. Кроме того, индонезийское золото родственно и сибирскому. Так что все методики, которые работают в Индонезии, равно применимы и в Мексике и в Сибири, хотя есть определенные нюансы из-за наличия слоя осадков в Сибири толщиной десятки метров, практически отсутствующего на рассматриваемом нами участке Западной Сумбавы. Впрочем, выходы рудных жил на поверхность не редкость и для Сибири, но такие выходы активно разыскивались и разрабатывались еще с царских времен, поэтому сейчас они чрезвычайно редки.

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

Теперь немного поговорим об использованном техническом подходе. Я уже не раз писал на хабре про построение модели плотности геологической среды по гравитационным данным (решение обратной задачи геофизики) и про уточнение глобальных гравитационных моделей низкого разрешения (для всей планеты) с помощью данных рельефа и космических снимков. В случае с наличием растительности оптические снимки мало пригодны, а вот радарные сигналы частично проходят сквозь растительность, особенно в сезон с минимумом растительности (или сухой растительностью). Собрав композит из множества разновременных радарных снимков, мы можем построить детальный рельеф (до 5 м по данным Sentinel-1 и даже детальнее, если снимков много), потом по рельефу детальную гравику и, наконец, геологическую модель высокого разрешения. Кстати, если использовать оптические снимки высокого разрешения для участков без растительности, можно и еще уточнять модель (PostgreSQL/PostGIS позволяют сделать такую обработку на разреженных данных, в моем гитхабе есть примеры).

Геологическое моделирование

Данные проб с концентрациями золота я взял из статьи, приведенной в списке литературы в конце статьи. Эта публикация интересна тем, что в ней приведены координаты набора скважин и детальные содержания по ним, хотя и без информации о глубинах. Впрочем, из картинки в статье очевидно, что речь идет о первых десятках метров. По причине высокой неравномерности концентраций золота вместо стандартного промышленного содержания золота в породе 3 г/тонну (3 ppm) будем считать значимым значение концентрации 5 г/тонну, так что у нас получится 4 продуктивных скважины и 12 не продуктивных. Заметим, что в дальнейшем, геолог полученные результаты проверит и убедится, что скважины с промышленным и полупромышленным содержанием приурочены к границам выделенных рудных зон.

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

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

Посмотрим на исследуемый участок вблизи:

Можно предположить, что скважины с высоким содержанием приурочены к естественным террасам. Отобразим модель геологической плотности отдельно:

Теперь стало понятнее золото накапливается в зонах низкой плотности, естественных карманах. Воспользуемся методом опорных векторов (SVC) с гауссовым ядром (помним про гауссовость модели геологической среды) для классификации наборов значений плотности по глубине для заданных классов продуктивности скважин. Чтобы "зацепить" структуры, по которым золотоносные растворы поднимались на поверхность, нужна модель глубиной минимум 300-500 метров (при меньшей глубине точность классификатора кардинально снижается по указанной причине). Найденные перспективные зоны оконтурены черной линией на поверхности:

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

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

Как мы помним, мы имеем дело с фрактальным строением, поэтому посмотрим более грубую модель с выделенными вертикальными аномалиями:

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

Заключение

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

Кстати, хотелось бы рассказать про золото в Сибири, там все интереснее, но даже карты СССР на эту территорию до сих пор засекречены. Самородное золото на поверхности практически выбрано (хотя до сих пор самородки находят), зато несколько метров вглубь и там начинается Эльдорадо. А серебра, палладия и платины и вовсе навалом, поскольку визуально они не выделяются. У меня уже припасено несколько интересных участков по соседству с теми местами, где я работал над проектами, так что надеюсь со временем кое-что показать.

Ссылки

• HOWTO: Multi-scale 3D Gravity and Magnetic Inversion modeling; https://www.linkedin.com/pulse/brief-explanation-3d-density-depth-model-construction-pechnikov/
• Gold grade of epithermal gold ore at Lamuntet, Brang Rea, West Sumbawa District, West Nusa Tenggara Province, Indonesia. AIP Conference Proceedings 1855, 050002 (2017); https://doi.org/10.1063/1.4985516 Rika Ernawati, Arifudin Idrus, and Himawan TBMP
• Sentinel-1 SAR GRD: C-band Synthetic Aperture Radar Ground Range Detected, log scaling; https://developers.google.com/earth-engine/datasets/catalog/COPERNICUS_S1_GRD
• Sentinel-2 MSI: MultiSpectral Instrument, Level-2A; https://developers.google.com/earth-engine/datasets/catalog/COPERNICUS_S2_SR
• ALOS DSM: Global 30m; https://developers.google.com/earth-engine/datasets/catalog/JAXA_ALOS_AW3D30_V3_2
• C-Support Vector Classification sklearn.svm.SVC; https://scikit-learn.org/stable/modules/generated/sklearn.svm.SVC.html

Эта статья является продолжением двух предыдущих: Ударим биспектром по бездорожью, или как найти золото в Сибири, в которой мы рассмотрели геологическую модель месторождения золота на территории Новосибирской области и Ищем рудное золото на острове Сумбава, Индонезия, в которых мы построили геологически обусловленную модель машинного обучения для поиска золота или других рудных минералов по всему Тихоокеанскому рудному поясу, используя для геологического моделирования открытые данные на платформе Google Earth Engine (GEE).

Убедившись в сходстве геологических моделей Сибири и Индонезии, применим классификатор для рудного золота Западной Сумбавы, Индонезия для Сибири. Так мы получим геологически корректный детальный прогноз золоторудности для Сибири, не используя вообще никаких геологических данных для этого региона.

Сопоставление моделей на территории Индонезии и Сибири

У нас есть уже обученный на детальных открытых данных классификатор для Индонезии, который мы хотели бы использовать и в Сибири. Сравним детальные модели плотности без учета влияния рельефа, построив для Западной Сибири модель с теми же параметрами, что и для Западной Сумбавы и выбрав равные участки. На левой модели точками на поверхности показаны скважины с высокими (Au) и низкими пробами золота (No Au), а на правой контурами отображены известные региональные месторождения золота (Au) и молибдена (Mo):

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

Классификация золоторудных территорий

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

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

Аналогично применим построенный классификатор к участку Западной Сибири (очевидно, здесь мы не можем оценить точность классификации автоматически):

Визуально результат соответствует ожидаемому и классификатор выделяет все те признаки, которые мы перечислили выше. Проверим на 3D модели:

Полученный результат геологически обоснован, поскольку выделенные области находятся в зонах накопления на границе выходов твердых пород, как мы это подробно разбирали ранее. При этом в области молибденового месторождения тоже выделены потенциальные золотые депозиты, что не удивительно, поскольку на региональной карте месторождения указаны по преобладающему минералу. Кроме того, на онлайн карте N-45-XIII. Схема прогноза полезных ископаемых, м-б 1:500 000 эти два месторождения объединены как золото-молибденовое месторождение.

Заключение

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

Если вас заинтересовало то, о чем я рассказываю, посмотрите ссылки ниже, вероятно, там вы найдете для себя интересное и полезное.

Ссылки

• Мои статьи на Хабре
• Теоретические и практические статьи и посты на LinkedIn
• Геологические модели и код на GitHub
• YouTube канал с геологическими моделями
• Геологические модели в виртуальной/дополненной реальности (VR/AR)

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

Видите взаимосвязь ортофотоснимка и рельефа? Если да, то вы или геолог или можете им стать: корреляция компонентов (разложения в пространственный спектр) составляет 41% для длины волны 20 м, 58% для 50 м и 99% для 300 м (Jupyter Python ноутбук с вычислениями доступен по ссылкам ниже). Большинство геофизиков поклянутся, что это у вас спектры порченые (записано с натуры), игнорируя и геофизику и прилагаемые вычисления и ссылки на публикации.

Миф первый, рельеф и сила тяжести

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

Важнейшим доказательством изостазии является отсутствие связи между рельефом и силой тяжести.

Увы, именно в отсутствие такой связи верят многие из встреченных мной геофизиков (но не геологи), хотя это абсолютно неверно! Что интересно, в англоязычной версии этой же викистраницы написано прямо противоположное, смотрите пример с айсбергом в разделе Deposition and erosion:

An analogy may be made with an iceberg, which always floats with a certain proportion of its mass below the surface of the water. If snow falls to the top of the iceberg, the iceberg will sink lower in the water. If a layer of ice melts off the top of the iceberg, the remaining iceberg will rise. Similarly, Earth's lithosphere "floats" in the asthenosphere.

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

Посмотрим, что говорит оглавление русскоязычного учебника по геофизике (Викулин, 2004):

Глобальные волны геоида, отсутствие их связи с особенностями строения земной коры

Здесь глобальные волны геоида это изменения (аномалии) поля силы тяжести порядка масштабов планеты, для которых, как и указано, отсутствует связь с рельефом. В то же время, для всех компонентов силы тяжести меньшего масштаба эта связь есть! Могу лишь предположить, что студенты-геофизики зачастую не понимают смысла выражения "глобальные волны геоида" (кстати, для обработки спутниковых данных современные глобальные модели оперируют "изменяющимся во времени квазигеоидом" с явным указанием порядка волны...).

Более того, корреляция компонентов силы тяжести (или вертикальной компоненты гравитационного поля) и рельефа близка к 100%, как показано в работе "Radially symmetric coherence between satellite gravity and multibeam bathymetry grids" (Marks and Smith, 2012). На графике ниже график радиальной когеренции (справа) из указанной статьи дополнен мною двумерной коррелограммой (слева) для глобальной модели гравитационных аномалий Sandwell and Smith Gravity Anomaly и глобальной модели рельефа GEBCO 2019:

Если на графике радиальной когеренции (корреляции) не понятно, почему для первых километров масштаба связь компонентов мала, то на двумерной коррелограмме все становится ясно измеренная сила тяжести относится к поверхности океана, а топография, очевидно, ко дну. Таким образом, здесь мы видим, что до масштабов ~40 км корреляция компонент близка к 100%. Более того, для региона Индонезия высокая корреляция (>75%) наблюдается до масштабов в сотни километров и средняя (>50%) для масштабов в тысячи километров, как показано в статьях по ссылкам ниже.

Миф второй, рельеф и аномалия Буге

Что интересно, миф второй полностью противоречит первому, что не мешает многим геофизикам верить в оба разом. Встречайте Аномалию Буге (Бугера):

Аномалия Бугера расхождение между теоретически вычисленным и реальным значением гравитационного поля Земли в определённой точке референц-эллипсоида.

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

Редукция Бугера называется простой или неполной, если местность можно приближенно представить в качестве бесконечной плоскости.

Очевидно, что с высотой (удалением от планеты) сила тяжести уменьшается, поэтому наблюдаемые (измеренные) на разной высоте значения силы тяжести нельзя сопоставить напрямую. Вот Буге и придумал способ, как удалить эту разницу из-за изменения высоты точки наблюдения (измерения), вычитая некую константу, домноженную на высоту точки наблюдения. Конечно, использование модели плоской Земли создает огромную ошибку при обработке современных точных данных, об этом я уже писал на Хабре в статье Почему в 21 веке геофизики верят в теорию плоской Земли?

И следующий шаг, или редукция:

Уточнённая, или полная редукция Бугера позволяет полностью учесть влияние рельефа местности.

Здесь речь о том, что наличие на местности неровностей в виде возвышений или углублений влияет на измеренную в этих местах силу тяжести, так что для сравнения измеренных значений нужно и их компенсировать. Заметим, что такая поправка, то есть необходимость учитывать влияние рельефа местности на силу тяжести, очевидно противоречит первому мифу. И эта поправка Буге проста в использовании нужно было среднюю плотность пород умножить на высоту точки наблюдения (возможно, Буге и до интеграла по глубине бы додумался, только где ему было взять среднюю плотность пород на территории в зависимости от глубины). Для самого Буге, с учетом расстояния между точками наблюдения километров 20 и более, о большой точности речь не шла, и средняя плотность собранных им же геологических образцов по региону его вполне устраивала для вычисления поправки. Увы, для современных детальных измерений этот метод совсем никуда не годится.

Посмотрим практический пример, как компоненты силы тяжести глобальной модели WGM2012 после редукции Буге (слева) и еще одной редукции (справа) коррелируют с компонентами рельефа глобальной модели GEBCO 2019:

Как видим, про "полностью учесть рельеф местности" речь тут явно не идет то есть результат редукции Буге противоречит его определению, а значит, его использование неправомерно. Обратим внимание на правый график тут мы ожидаем максимум корреляции 100%, а наблюдаем всего около 75%, что показывает неточность самой модели WGM2012 (по сравнению с использованной выше Sandwell and Smith Gravity Anomaly, авторы которой вовсе не вычисляют редукцию Буге по уже понятным причинам), таким образом, здесь реальная корреляция редукции Буге с рельефом получится около -100%.

Если обратиться к учебнику "Гравиразведка" (Утёмов, 2009), там в соответствующем разделе есть разъяснение полученного результата:

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

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

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

Миф третий, сферический конь эллипсоид в вакууме

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

В геофизике модель геопотенциала представляет собой теоретический анализ измерения и расчета эффекты Земли гравитационного поля

Да-а, расчета эффекты Земли это, конечно, не совсем русский язык, или даже совсем не русский. Равно как и теоретический анализ измерения тоже тот еще шедевр. Тем не менее, можно понять, что модель геопотенциала это модель (глобального) гравитационного поля Земли. Это определение уже ошибочное на самом деле, рассматривается вертикальная компонента гравитационного поля, то есть сила тяжести (гравитационное поле векторное, а сила тяжести скалярная, и связь между ними весьма сложная). Но, главное, геофизикам такая модель вообще не нужна! А кому нужна, описано дальше в этой же викистатье:

Рекурсивные алгоритмы, используемые для численного распространения орбит космических аппаратов

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

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

Миф четвертый, о пространственных спектрах

У меня в коллекции есть такие перлы, как спектры у вас порченые (освятить, что ли, хотят?) и другие. Попросту говоря, геофизики не имеют никакого понятия, что это такое.

Давайте посмотрим, как же практически работать с пространственными спектрами. Как всегда в математике, все эти преобразования по сути чрезвычайно просты. Например, результат применения гауссова фильтра масштаба N метров есть низкочастотная фильтрация (low-pass filter), разница между исходным изображением и фильтрованным есть высокочастотная фильтрация (high-pass filter) и разница между двумя фильтрами масштаба N и M метров есть полосовая фильтрация (band-pass filter). Для выделения компонента пространственного спектра шириной l метров для длины волны L метров достаточно выполнить полосовую фильтрацию с фильтрами масштаба L-l/2, L+l/2. Спектр мощности может быть представлен как стандартное отклонение полученных компонент, а через логарифм спектра мощности посчитана и фрактальная размерность. Для примера посмотрите график радиальной когерентности из статьи НАСА и двумерную спектрограмму, полученную указанным методом, в обсуждении первого мифа выше.

Низкочастотная фильтрация это улучшенная замена Буге преобразования:

А если выполнить высокочастотную фильтрацию, то мы сразу получаем все глобальные разломы то есть границы литосферных плит и микроплит:

У специалистов может возникнуть вопрос, зачем стоит использовать (относительно медленную) матричную свертку для вычислений вместо, скажем, быстрого Фурье-преобразования (БПФ). Так вот, плюсы такого подхода заключаются в том, что при вычислении свертки мы можем корректно обрабатывать пропущенные значения (игнорировать их, а не использовать значение-заполнитель, и одновременно подсчитывать количество, считая результат неопределенным при большом числе пропущенных значений) и границы и даже работать с разреженными данными.

Интересно, что аргентинские геофизики применяют индекс фрактальности для вычисления средней плотности по территории в зависимости от глубины для использования в преобразовании Буге. В самом деле, этот способ работает для территории Латинской Америки (и не только там, но не везде по миру все зависит от характера профиля плотности) я проверил с помощью численного моделирования, по ссылкам ниже вы найдете готовый Jupyter Python ноутбук с вычислениями. А вот глобальные модели гравики для Латинской Америки с помощью сферических разложений получаются очень так себе (впрочем, для космических вычислений это не мешает, там не нужны детальные пространственные компоненты силы тяжести на уровне поверхности) слишком высока фрактальность территории (сильно расчлененные горы с их отрогами), видимо, потому и пришлось придумывать фрактальный метод анализа.

Заключение

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

Бонус

И специальное дополнение для мифоборцев. Не стал добавлять еще один миф, просто "на закуску" приведу выдержку из рецензии, когда-то полученной от журнала "Отечественная геология":

Авторы пытаются нас убедить в том, что математические манипуляции с оптическим сигналом одного из каналов мультиспектральной съемки позволяют нам анализировать глубинное строение Земли. Математические обработки мультиспектральных снимков направлены на улучшение их визуального восприятия и классифицирования представленных на них объектов для создания тематических карт.

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

Ссылки

Marks, K.M., Smith, W.H.F., 2012. Radially symmetric coherence between satellite gravity and multibeam bathymetry grids. Mar Geophys Res 33, 223227. https://doi.org/10.1007/s11001-012-9157-1

Викулин А.В., 2004. ВВЕДЕНИЕ В ФИЗИКУ ЗЕМЛИ. УЧЕБНОЕ ПОСОБИЕ ДЛЯ ГЕОФИЗИЧЕСКИХ СПЕЦИАЛЬНОСТЕЙ ВУЗОВ

Утёмов Э.В., 2009. Гравиразведка

Spectral Coherence between Gravity and Bathymetry Grids

WGM2012 spatial components of free-air gravity and topography correlation

Bouguer and Free-Air Gravity Anomalies in terms of spatial spectrum

There is a high correlation between DEM and ortho photos

3D Density Inversion by Circular Hough Transform [Focal Average] and Fractality Index

[Gaussian Filtering on Spheroid [Sandwell & Smith]](

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

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

Измерение фрактальной размерности

Существует множество способов определения фрактальной размерности, и общего у них только наличие характерного пространственного масштаба. Можно производить вычисления фрактальности в пространственной области, а можно и в частотной (компоненты пространственного спектра), дифференцировать и интегрировать Секрет в том, что всё многообразие способов дает схожие результаты и выбор конкретного способа довольно произволен, важно лишь понимать, как результаты выбранного метода вычислений сопоставить с результатами других методов (Gneiting et al., 2012).

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

Производная (или дифференциал) означает наклон или приращение кривой, а интеграл (или отрицательная производная) означает заполнение или площадь под кривой. Таким образом, для произвольной кривой (одномерного сигнала) дифференциальный способ оперирует ломаной, характеризующей границу кривой, а интегральный способ площадью под кривой. Удивительно здесь то, что и дифференцирование и интегрирование задают пространственный масштаб, или окно вычислений. С интегралом сразу понятно, ведь площадь, очевидно, зависит от ширины фигуры под кривой, а вот с дифференциалом все не так просто. Подумаем, сколько значений минимально необходимо для вычисления первого дифференциала. Очевидно, два ведь нам нужно посчитать разницу между ними. А для второго дифференциала нам нужно вычислить разницу между двумя первыми дифференциалами, каждый из которых требует два значения и, если взять одно общее значение для вычисления первых дифференциалов, то потребуется три значения. Как видим, для вычисления N-го дифференциала необходимо N+1 значений. Значит, степень производной (дифференциала) также определяет пространственный масштаб. Таким образом, можно определить фрактальную размерность через интегрирование по отрезку изменяемого размера (окну), или через производные разных порядков (включая дробные и отрицательные, при этом последние соответствуют одной из первообразных, то есть интегралам). За подробностями рекомендую обратиться к публикациям ныне покойного профессора Нижегородского государственного университета имени Н. И. Лобачевского (Университет Лобачевского, ННГУ) Александра Ивановича Саичева, в которых я и нашел связь фрактальности с дробными производными (Саичев, В. А. Филимонов, 2008).

Корреляционная размерность и (поли)спектры определяются наборами значений с заданными расстояниями между ними, то есть содержат пространственный масштаб. Через них, соответственно, тоже можно определять фрактальную размерность. Основой для описания и понимания (поли)спектрального и корреляционного видов анализа является анализ кумулянтный (Дубков, Малахов, 1973). Для детального ознакомления рекомендую к прочтению теоретические и практические работы по кумулянтному и биспектральному (полиспектральному) анализу от преподавателей ННГУ Александра Александровича Дубкова и Германа Николаевича Бочкова, у которых мне посчастливилось учиться. Пользуюсь случаем поблагодарить их обоих, а Германа Николаевича еще и поздравить с юбилеем 80 лет!

Далее мы будем пользоваться вычислением фрактальной размерности через дисперсию компонентов оконного дискретного спектра мощности. В пространственной области аналогичный анализ выполняется круговым преобразованием Радона с дисперсией как базовой функцией. Иными словами, если в ранее нами использованном круговом преобразовании Радона, см. Методы компьютерного зрения для решения обратной задачи геофизики, заменить вычисление среднего значения на дисперсию (стандартное отклонение), тоже получится способ вычисления фрактальной размерности, притом довольно популярный (хотя и под разными названиями) в геологической литературе (см. ссылки в конце статьи). В радиофизике для описания подобных преобразований есть так называемые методы синтеза апертуры, то есть получения большого количества отсчетов вместо одного. Классический пример локатор бокового обзора (эхолот, в том числе) излучает и принимает отраженные сигналы непрерывно и для вычисления расстояния до каждой точки береговой линии или рельефа дна используется множество измерений, что позволяет сильно улучшить точность. Для нашей задачи рассмотрим простой пример можем ли мы определить фрактальность в одной точке рельефа разрешением 10 м? Очевидно, нет. Если же мы проанализируем все точки в кольцах радиуса от 1 до 1000 пикселов (в диапазоне 10 м 10 км для заданного разрешения), включающих от 4х до ~6000 пикселов (2R, где R это радиус кольца), то сможем с высокой точностью вычислить фрактальную размерность. При этом точность вычисления пропорциональна квадратному корню из числа пикселов или радиуса колец. А теперь в игру вступают силы природы если мы наблюдаем фрактальность на масштабе от 1 пиксела и более, эта же фрактальность гарантированно присутствует и на меньших масштабах, хотя мы не можем определить точно граничный масштаб. Небольшое отступление: в лабораторных экспериментах с моделированием динамики фотополимеров мы находили этот граничный масштаб в области сотен нанометров (где молекулы образуют сетчатые структуры и не гауссовые и не фрактальные), так что для геологических исследований граничного масштаба просто нет. Таким образом, благодаря самой сути явления фрактальности, становится возможным обнаружить области залегания фрактальных рудных тел размером значительно менее метра при анализе спутниковых оптических или радарных снимков 5-10 м разрешения (например, Sentinel-1 и Sentinel-2). Смотрите практические примеры в публикации Ударим биспектром по бездорожью, или как найти золото в Сибири. Обратите внимание, что сегодня мы рассматриваем спектральный метод, а не биспектральный, который обладает дополнительным важным свойством исключения гауссовых сигналов и, таким образом, значительно более чувствителен за счет полного игнорирования всех гауссовых помех, представляющих неразрешимую проблему для спектрального анализа на субпиксельных масштабах.

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

Изостазия и граница фрактальности рельефа и силы тяжести

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

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

ИЗОСТАЗИЯ (от изо и греч. состояние), состояние механич. равновесия поверхностной оболочки Земли, при котором на некоторой глубине (обычно принимаемой равной 100150 км) в верхней мантии давление вышележащих пород становится одинаковым. В этом случае избыток или недостаток масс на поверхности Земли компенсируется соответствующим перераспределением масс в её недрах. [Большая Российская Энциклопедия Изостазия]

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

кора плавает в вязкой мантии в соответствии с законом Архимеда, т. е. находится в состоянии гидростатич. равновесия. [Большая Российская Энциклопедия Изостазия]

С учетом известного нам соотношения между длиной волны и глубиной ее гравитационного источника, см. Методы компьютерного зрения для решения обратной задачи геофизики, для пространственных компонент рельефа и поля силы тяжести верхняя граница фрактальности составляет примерно 150-300 км, варьируясь для разных регионов планеты. Вполне очевидно, что (сжатая вязкая) жидкость не обладает свойством фрактальности, в отличие от твердых тел, что и является причиной исчезновения эффекта.

Что интересно, известно об эффекте изостазии уже почти три столетия, и первооткрывателем был тот самый Буге (Бугер), результаты работ которого и сегодня вызывают столько недопонимания, см. Легенды и мифы геофизики:

Данные о том, что вес горного сооружения компенсируется более лёгкими массами на глубине, впервые получил П. Бугер в 1749. Проанализировав результаты экспедиции 1736, он обнаружил, что в предгорьях Анд угол уклонения отвеса от вертикали значительно меньше того, который должны создавать массы горного рельефа. [Большая Российская Энциклопедия Изостазия]

Определение границы Мохоровича (Мохо) по значениям фрактальности и корреляции пространственных спектров рельефа и силы тяжести

Обратимся снова к вышеупомянутой статье энциклопедии:

Амер. геолог Дж. Баррелл в 191415 и нидерл. геофизик Ф. Вейнинг-Мейнец в 1931 развили эту схему изостатич. компенсации. Они предположили, что верхняя часть литосферы является упругой, и использовали решение задачи об изгибе внешними силами упругой плиты, плавающей на жидком основании. [Большая Российская Энциклопедия Изостазия]

Из вышесказанного становится понятно, что фрактальность проявляется максимально для твердого слоя (земной коры), далее резко снижается для слоя упругого и совсем исчезает при переходе к вязкой магме. Граница, или подошва, твердой земной коры определяется так называемой границей Мохоровича (Мохо), и до этой границы наблюдаются высокая фрактальность и корреляция пространственных спектров силы поля тяжести и рельефа близка к единице (100%).

Посмотрим на графики из публикации Легенды и мифы геофизики:

Здесь правый график является диагональным сечением левого. Смещение в начале оси абсцисс на графиках соответствует глубине океана на рассматриваемой территории (2.5 5 км). Очень высокая корреляция (90% и выше) наблюдается до маштаба ~15 км, высокая корреляция (75% и выше) до ~35 км, и дальше наблюдается плавный спад значимости корреляции. Как мы рассмотрели ранее, масштаб, на котором исчезает значимая корреляция, и будет соответствовать достижению изостазии. Для океанского дна на рассматриваемой нами территории граница Мохо расположена на глубине 11 км согласно модели CRUST 1.0: A New Global Crustal Model at 1x1 Degrees, что соответствует длине волны 15 км, как мы и видим на графиках. Граница слоя упругости составит около 25 км (длина волны 35 км). Таким образом, по корреляции компонентов пространственных спектров поля силы тяжести и рельефа можно определить границы Мохоровича (Мохо), упругого слоя и изостазии.

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

Горизонтальное смещение кривых плотности вызвано приведением поля силы тяжести к поверхности, а батиметрии, очевидно, ко дну океана. Здесь мы также легко выделяем области пространственных масштабов около 15 км и 35 км, как и на предыдущих графиках. Как будет показано далее, значение фрактальной плотности по данным рельефа является оценкой сверху для геологической плотности и при увеличении площади наблюдения при ее неоднородности приводит к систематическому завышению оценки плотности фрактальным методом.

Сила Кориолиса и фрактальность рельефа и силы тяжести

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

Если мы рассмотрим острова Индонезии в Южном полушарии, то потоки гидротермальных растворы с минералами под действием силы Кориолиса должны отклоняться влево от направления движения. Посмотрим на региональную модель из статьи Ищем рудное золото на острове Сумбава, Индонезия:

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

В Западной Сибири, расположенной в Северном полушарии, гидротермальные потоки отклоняются вправо. Вернемся к модели из статье Ударим биспектром по бездорожью, или как найти золото в Сибири:

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

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

И снова о редукции Бугера (Буге)

В уже упомянутой выше публикации Легенды и мифы геофизики мы показали, что на самом деле, редукция Буге не удаляет (полностью) корреляцию спектральных компонентов измеренной силы тяжести и рельефа. Стоит отметить, что основная проблема проявляется для диапазона масштабов до 20 км, что вовсе не было проблемой для самого Буге, поскольку детальность его данных была ниже и для них придуманная им редукция работала практически идеально. Поскольку это просто классический секрет Полишинеля, то придуманы и способы исправить ситуацию например, исходя из равной фрактальности рельефа и силы тяжести, можно подобрать параметры преобразования Буге, чтобы несколько уменьшить фрактальность (Miranda et al., 2015). В работе (Zhang, Featherstone, 2019) для территории Австралии методами фрактального и спектрального анализа показано, что топография, редукция Буге и аномалии в свободном воздухе имеют почти идентичные фрактальные размерности и спектры, при этом использование редукции Буге приводит к наибольшим ошибкам.

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

Фрактальность рельефа, силы тяжести и снимков

При моделировании поля силы тяжести от набора гауссовых глубинных источников 3D Density Inversion by Circular Hough Transform (Focal Average) and Fractality Index оказывается, что это поле на поверхности обладает выраженной фрактальной размерностью, причем эта размерность пропорциональна распределению плотности с глубиной. Таким образом, становится ясно, что первопричиной является фрактальность гравитационного поля, которая с высочайшей корреляцией проявляется и в рельефе, следовательно, эрозионные процессы не нарушают эту взаимосвязь. Для космических или ортофотоснимков значение фрактальности оказывается выше за счет вариаций цвета, не связанных с изменениями рельефа. Выбирая подходящий спектральный диапазон, в котором поверхность имеет наименьшую спектральную размерность, мы можем оценить плотность пород поверхности. Обратимся к доступному на GitHub ноутбуку Geological Fractality on ALOS DEM AW3D30 and Sentinel-2 SUrface Reflectance Image for Ravar, Kerman Province, Iran:

Плотность, вычисленная по фрактальной размерности рельефа, равна =3200 kg/m при R=0.99, что вполне согласуется с плотностью магматических пород, образующих эту территорию. По космическому снимку значение плотности выше и составляет =3500 kg/m при R=1.00, что соответствует плотности поверхностных метаморфических пород, вероятно покрывающих часть территории. При этом следует учитывать, что индекс фрактальной размерности рассчитывается по вариации значений и систематически завышает оценку в том случае, если территория неоднородна по составу. Таким образом, полученные нами значения являются оценкой сверху. Для получения более точных значений необходимо сегментировать территорию на гомогенные участки и выполнить вычисления для них по отдельности. Кроме того, перегибы кривой фрактальной размерности пропорциональны глубинам, на которых геологическая плотность изменяется.

Заключение

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

Также смотрите

• Мои статьи на Хабре
• Теоретические и практические статьи и посты на LinkedIn
• Геологические модели и код на GitHub
• YouTube канал с геологическими моделями
• Геологические модели в виртуальной/дополненной реальности (VR/AR)

Ссылки

Кумулянтный анализ функционального нелинейного преобразования негауссовых случайных процессов и полей, А. А. Дубков, А. Н. Малахов, 1973. http://m.mathnet.ru/links/cd623304046883b7c36697e2e9f9b1d0/dan39067.pdf

Кумулянтный анализ случайных негауссовых процессов и их преобразований, Малахов А.Н.,1978. https://ikfia.ysn.ru/wp-content/uploads/2018/01/Malahov1978ru.pdf

Численное моделирование реализаций и спектры квазимультифрактального диффузионного процесса, А. И. Саичев, В. А. Филимонов, 2008. http://www.mathnet.ru/links/2a351947644994381a8272c4fc3ba0dd/jetpl108.pdf

Numerical simulation of the realizations and spectra of a quasi-multifractal diffusion process, A. I. Saichev & V. A. Filimonov, https://link.springer.com/article/10.1134/S0021364008090129

Gneiting, T., evkov, H. & Percival, D. B. Estimators of Fractal Dimension: Assessing the Roughness of Time Series and Spatial Data. Statist. Sci. 27, 247277 (2012).

Zhang, K. & Featherstone, W. Exploring the Detailed Structure of the Local Earths Gravity Field Using Fractal and Fourier Power Spectrum Techniques. (2019).

Miranda, S. A. et al. Fractalness of land gravity data and residual isostatic anomalies map of Argentina, Chile and western Uruguay. Geofsica internacional 54, 315322 (2015).

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

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

Спутниковая карта Google Satellite наложена на детальный рельеф USGD NED DEM 1 м

Введение

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

Рассмотрим ниже некоторые принципиальные особенности цифрового рельефа и далее проверим его качество на примерах от различных производителей.

Заявленное пространственное разрешение рельефа

Пространственное разрешение рельефа, даже построенного по одним и тем же исходным данным, варьируется в широких пределах в зависимости от территории. И, тем более, разное разрешение и точность имеют результаты, полученные разными способами скажем, методом триангуляции космоснимков и с помощью лидарной съемки. К примеру, для территории США доступно множество продуктов рельефа разрешением от 1 м до 30 м, так что вся территория страны покрыта рельефом 30 м и 10 м, а часть территорий доступны с разрешением 5 м и 1 м. Таким образом, единый рельеф разрешением 1 м на всю территорию США будет синтезом разномасштабных рельефов и получившиеся детализация и точность будут варьироваться по территории. В идеале, следовало бы объединять данные в спектральной области или с использованием интерполяции методом ближайшего соседа, на практике же часто используются нелинейные методы интерполяции, так что получившийся продукт содержит широкую полосу мусорных компонент пространственного спектра. При взгляде на такой рельеф становится понятно, что выглядит он как-то не так, но точную оценку можно получить лишь при анализе его пространственного спектра.

Заявленная вертикальная точность рельефа

Точность рельефа может определяться совершенно разными способами, например, как величина ошибки относительно набора референсных точек на поверхности или относительно исходных данных (вопрос точности которых это совсем другая история). Точность может указываться и так, к примеру: ошибка не более 5 м, что на самом деле означает ошибку не более 5 м с доверительным интервалом 95% (или другим), то есть вовсе не гарантирует точность 5 м для любого отдельно взятого пиксела или участка. Поскольку точность оценивается для отдельных пикселов и отдельных участков, для которых есть точные отметки высот, то в пределах большой территории может сильно варьироваться. Например, если 99% рельефа занимает плоская равнина с малыми перепадами высот и, следовательно, высокоточным рельефом, то оставшийся 1% рельефа может иметь точность в 100 раз худшую. Поэтому рельеф заявленной точности 5 м доступен на всю территорию планеты, а 10 см точности только выборочно. Но и это еще не все. Точность рельефа видимой поверхности (Digital Surface Model, DSM) соотносится к точности использованных для создания рельефа данных к измеренной поверхности (в зависимости от местности это может быть лес или скалы и так далее), так что при другом ветре и в другой сезон измеренные значения окажутся далеко за пределами заявленной точности. В случае же рельефа непосредственно поверхности планеты (Digital Terrain Model, DTM) есть разные методы исключения растительности (даже трава и кустарник дают погрешность по высоте более 10 см, не говоря про деревья), а оценка точности производится, как правило, по некоторым референсным точкам только на открытой местности.

Оценка реального пространственного разрешения рельефа

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

Оценка вертикальной точности рельефа

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

Глобальный рельеф всей планеты ALOS World 3D 30m (AW3D30) вертикальной точностью 5 м, построенный методом триангуляции. Реальное разрешение 30 м

Это комбинированный открытый продукт размером около 220 ГБ, доступный на сайте производителя ALOS World 3D 30m (AW3D30) и на платформе Google Earth Engine (GEE) как ALOS DSM: Global 30m. Комбинированный он потому, что использует для заполнения пропущенных значений рельефы SRTM 30m, ASTER DEM и другие. На мой взгляд, является лучшим из открытых глобально доступных. Если SRTM содержит серьезные пиксельные артефакты, а ASTER DEM буквально кляксы некорректно интерполированных значений, то ALOS практически не грешит подобными проблемами. Анализ пространственного спектра в Python ноутбуке показывает следующий результат:

Для пространственного спектра рассматриваемого рельефа в двойных логарифмических координатах коэффициент детерминации R-квадрат равен 98% для масштаба от 30 метров. Таким образом, этот рельеф явно получен из более детального и действительно, поставщик также предлагает коммерческий рельеф разрешением 5 м и все той же вертикальной точностью 5 м.

Коммерческий рельеф 1 м, построенный методом триангуляции. Реальное разрешение skipped

Образец коммерческого рельефа со спутникового аппарата Pliades от компании Airbus DS Geo SA получен с официального сайта Elevation1 DSM + Pliades Ortho 0.5m pan-sharpened (Orthoimage included). Для целей тестирования я выбрал участок горной местности почти без техногенных объектов. Обратим внимание на лицензионное соглашение, которое разрешает только внутреннюю техническую оценку продукта ("to use the PRODUCT for internal technical evaluation purposes only") и запрещает публикацию любых результатов ("any derivative product or information"). В связи с этим, я опубликую только Python ноутбук, который вы можете запустить для собственной "внутренней технической оценки продукта", согласно лицензии, и получить точное значение пространственного разрешения. Обратите внимание, что отметка 60 м на графике пространственного разрешения в коде ноутбука поставлена исключительно для удобства оценки фрактальности спектра и я не несу ответственности за то, что она может вам показаться равной реальному пространственному разрешению рассматриваемого рельефа.

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

Рельеф США USGD NED DEM 1 м вертикальным разрешением 10 см, построенный по данным лидарной съемки. Реальное разрешение около 2 м

Предоставляется на Amazon AWS Staged Products Directory в виде ~7700 исходных zip архивов общим объемом около 1 ТБ и распакованным размером 1.7 ТБ. Реальная карта покрытия составляет примерно половину от официально заявленной я скачал все тайлы и сделал карту покрытия, на которой оказалось, что многие тайлы целиком или существенно перекрываются (хотя не должны, по документации перекрытие соседних тайлов составляет ровно 6 граничных пикселов). Содержит заметные пиксельные артефакты, но качество впечатляет по сравнению с рельефом, построенным методом триангуляции снимков. Анализ пространственного спектра в Python ноутбуке показывает следующий результат:

Для пространственного спектра рассматриваемого рельефа в двойных логарифмических координатах коэффициент детерминации R-квадрат равен 98% для масштаба от 2х метров (для масштаба от 4х метров R-квадрат равен 99%). Таким образом, как и следует из теории, точный рельеф фрактален на всех масштабах.

Заключение

С помощью качественного рельефа высокого разрешения можно выполнить гидродинамическое моделирование на поверхности, смотрите Гидродинамическое моделирование (CFD) на рельефе с помощью MantaFlow и визуализация результатов в ParaView, и построить детальные геологические модели как показано в статье Построение достоверных геологических моделей. При использовании методов фрактальной математики становится возможным выделение рудных объектов метрового масштаба, см. Пространственные спектры и фрактальность рельефа, силы тяжести и снимков. Также детальный рельеф помогает, при определенных условиях, уточнить спутниковые интерферограммы и получить более детальную картину отражения сейсмических волн от глубинных объектов и более точную модель смещений поверхности, см. Геология XXI века как наука данных о Земле и Вычислительная геология и визуализация.

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

Также смотрите

• Мои статьи на Хабре
• Теоретические и практические статьи и посты на LinkedIn
• Геологические модели и код на GitHub
• YouTube канал с геологическими моделями
• Геологические модели в виртуальной/дополненной реальности (VR/AR)

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

Чтобы начать редактирование в JOSM, необходимо скачать jar-файл с сайта josm.openstreetmap.de и установить его. В окне настроек можно указать язык, настроить авторизацию и скачать необходимые модули. Одними из наиболее удобных и широко используемых модулей являются:

• Buildings_tools инструмент для прорисовки зданий. Позволяет рисовать прямоугольники или круги, нарисованной фигуре проставляется по умолчанию тег building=yes;
• Utilsplugin2 добавляет несколько утилит, упрощающих редактирование геометрических фигур, а также расширяющих возможности выборки;
• Reltoolbox упрощает процесс работы с отношениями;
• Imagery_offset_db позволяет получить/сохранить ближайшее к месту редактирования смещение снимка;
• HouseNumberTaggingTool позволяет достаточно быстро тегировать адресную информацию для зданий, номера которых идут с определенным шагом.

Подробную информацию о настройке и начале работы с JOSM можно почерпнуть из видео на сайте josm.ru.

Итак, мы сделали первоначальную настройку редактора. Теперь подготовим рабочую область для рисования. Сперва надо в левом верхнем углу интерфейса найти кнопку с зеленой стрелкой вниз она позволит выбрать область для редактирования. В открывшемся модальном окне ставим галочки в разделе Источники и типы данных, чтобы скачать данные OpenStreetMap, GPS-треки и заметки. Прямоугольником на карте выделяем область, которую хотим редактировать (перемещение по карте в JOSM осуществляется с зажатой правой клавишей мыши) и кликаем на кнопку Скачать.


Если была выделена слишком большая область, JOSM выдаст сообщение о неудачном запросе (но GPS-треки и заметки будут загружены).


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


В правой части экрана располагается панель слоёв, в которой можно управлять их видимостью и порядком отображения. Добавим слой со спутниковыми изображениями: в главном меню развернем пункт Слои и выберем слой, например, Снимки Maxar Premium (Beta). На экране и в списке слоев появится слой со спутниковыми изображениями.


Настроим смещение снимка, кликнув в верхней панели инструментов на иконку с двумя красными стрелками.


Теперь рабочая область готова к рисованию. Начнем, как и в редакторе iD, с добавления здания. В JOSM есть два основных режима редактирования: режим выбора горячая клавиша S или верхняя кнопка Выделять, перемещать, масштабировать и вращать объекты на панели инструментов слева; и режим рисования горячая клавиша A или кнопка Рисовать точки на панели инструментов слева.

Здание часто представляет собой замкнутый многоугольник. В большинстве случаев жилой дом имеет форму, близкую к прямоугольной. Чтобы нарисовать ровный прямоугольник при установленном модуле Building_tools, достаточно нажать горячую клавишу B (от слова Building) и щелчками мыши на экране растянуть длину и ширину прямоугольника. Если же этот плагин не установлен, то можно начать редактирование горячей клавишей A (от слова Add), нарисовать многоугольник (двойной щелчок заканчивает рисование) и нажать горячую клавишу Q, которая спрямит углы. То же самое можно сделать, не пользуясь горячими клавишами, а нажимая на соответствующие кнопки в интерфейсе.


Более сложные формы можно рисовать несколькими прямоугольниками, выделить их все, сменив клавишей S режим редактирования на режим выделения, и горячими клавишами Shift+J объединить в один полигон.


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


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


Основной тег для многоугольника здания building. В соответствии с назначением здания необходимо проставить корректное значение этого тега. В простых случаях назначение здания совпадает с его использованием. Так, для многоквартирных домов используется значение apartments, для частных жилых домов detached, школы обозначаются school, а детские сады kindergarten. Чуть более замысловатым является тегирование зданий, использование которых не совпадает с их назначением: пример, рассмотренный в предыдущей статье здание, построенное как склад, в настоящее время используется как жилой дом. В этом случае здание тегируется как склад (building=warehouse), текущее использование задается тегом building:use=residential. Как быть, если, например, бывшее здание больницы используется теперь как библиотека? Ответ прост: здание должно иметь теги building=hospital и amenity=library. В сложных случаях, когда назначение дома определить не удается, применяется тегирование building=yes.

Итак, тип объекта был установлен. Теперь нужно указать адресную информацию. Она обозначается тегами с префиксом addr:. Минимально необходимая адресная информация для дома это его номер и название улицы, обозначающиеся тегами addr:housenumber и addr:street соответственно. При указании номера дома, состоящего не только из цифр, следует руководствоваться следующими правилами: литеры пишутся слитно с номером дома, а номера корпусов, строений, сооружений, владений и т.п. указываются через пробел. Например, номер дома 8А, корпус 1, строение 3 должен быть указан как 8А к1 с3.


Наиболее распространенным типом адресации является адресация по улицам, но как быть в случае, когда адрес относится к определенной территории, как, например, в Ангарске, Набережных Челнах или Троицке? В таких случаях необходимо использовать тег addr:place, который заменяет тег addr:street (то есть недопустимо одновременное использование тегов addr:street и addr:place для одного объекта).

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

Хорошим тоном будет указать почтовый индекс, задав тег addr:postcode.

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


На рисунке выше видно, что дома 14 и 18 имеют два подъезда, а в доме 16А их не менее пяти.

Вход в подъезд обозначается как entrance=staircase (за таким входом обычно скрывается лестница), черный вход entrance=service, пожарный выход обозначается как entrance=emergency, вход в частный дом или квартиру с улицы entrance=home. Если есть сомнения, как обозначать конкретный вход, допускается использование общего значения тега yes. Порядковый номер входа указывается с помощью тега ref. Для указания номеров квартир используется тег addr:flats. Указание диапазона квартир производится с помощью дефиса: addr:flats=1-28.



В нестандартных ситуациях, когда нумерация прерывается, диапазоны квартир в подъезде указываются через точку с запятой: addr:flats=1-12;14;16-20.

Дополнительно можно указать инженерную информацию о здании этажность (тег building:levels), форму крыши (roof:shape), материал фасадов здания и его крыши (теги building:material и roof:material соответственно), цвет внешних стен и крыши здания (теги building:colour и roof:colour). Однако новичкам следует с осторожностью указывать такую информацию для сложносоставных зданий, имеющих различную этажность и меняющуюся форму крыши. Для таких зданий в первом приближении можно отрисовать общий контур здания, проставить максимальную этажность. Для внесения более точных характеристик здания необходимо разделять его на элементарные части building:part, из которых собирается финальное отношение c типом multipolygon с тегом building. Вот, например, как выглядит собор Василия Блаженного в Москве:


Возможно, кому-то покажется, что отношения это очень сложно. Это и впрямь непросто, но стоит их освоить, и станет понятно, насколько они удобны и упрощают жизнь картографам. Разберем один случай, не такой изощренный, как представленный выше собор. Для работы понадобятся модули utilsplugin2 и reltoolbox. Удостоверьтесь, что они установлены. Наш случай это многоквартирный дом переменной этажности в жилом квартале Влюблино в Москве. Находится он по адресу улица Цимлянская, дом 3 корпус 2. Дом состоит из двух отдельно стоящих зданий, которые имеют секции с разной этажностью. Начнем с того, что нарисуем контуры отдельных зданий.


Теперь прорисуем линии, разделяющие разноэтажные секции дома. Выделив полигон здания и входящие в него перемычки и нажав горячую клавишу Q, можно скорректировать прямые углы, если это необходимо. Следующий шаг разбить цельные полигоны теми перемычками, которые мы нарисовали. Для этого выделяем полигон и входящую в него перемычку (важно делать это по одной) и нажимаем горячие клавиши Alt+X, либо на вкладке Еще инструменты находим пункт Разрезать объект.


Выделим получившиеся 5 полигонов и проставим им тег building:part=yes. Добавим также тег, описывающий форму крыши roof:shape=flat. Теперь, выделяя каждый отдельный полигон, проставим ему количество этажей building:levels, используя фотографические изображения зданий. Для трехсекционного здания это будут 1, 9 и 17; для двухсекционного 9 и 17.



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

Выделим снова все пять полигонов. В панели Relation Toolbox, расположенной в правой части экрана, нажмем на кнопку Multi. Эта операция разбила полигоны на отдельные сегменты и создала отношения на месте бывших полигонов building:part. При этом все заданные теги сохранились в соответствующих отношениях. Соберем итоговый мультиполигон для всего здания, он должен соответствовать тем контурам, которые мы нарисовали в самом начале. И поскольку два отдельных здания относятся к одному адресу, они должны быть участниками одного отношения (в нашем случае одного мультиполигона). Выделив линии, образующие два внешних контура здания, снова нажимаем кнопку Multi в панели Relation Toolbox и видим появившуюся у контуров зданий желтую обводку.



В панели Теги добавляем к уже имеющемуся тегу type=multipolygon теги building=apartments, building:levels=17, addr:street=Цимлянская улица, addr:housenumber=3 к2.

Описанный способ не является единственно верным, возможности OpenStreetMap позволяют один и тот же случай отрисовать несколькими способами. Если нет уверенности, как правильно нарисовать сложное многоэтажное здание, лучше задать вопрос в Telegram-сообществе OpenStreetMap, где бывалые осмеры подскажут, как быть в данной ситуации, и поделятся ссылками на места в интернете, где описаны правила или договоренности на интересующую тему.

Итак, первый этап пройден, дом нарисован. Теперь хорошo бы нарисовать окружающую инфраструктуру. Проезды, обозначенные знаками Жилая зона, имеющие название, тегируются как highway=residential + living_street=yes, неименованные внутридворовые проезды обозначаются как highway=service + living_street= yes. Пешеходные дорожки тегируются highway=footway. Детские площадки обозначаются leisure=playground, площадки для выгула собак leisure=dog_park, спортивные площадки leisure=pitch, желательно с указанием вида спорта (тег sport), парковки amenity=parking.

Когда с оцифровкой закончено, наступает время отправить свои правки на сервер OpenStreetMap. Для этого на верхней панели инструментов находим кнопку с зеленой стрелкой вверх. В этот момент происходит валидация сделанных правок, и если что-то сделано некорректно, JOSM выведет об этом сообщение.


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


В новом окне необходимо указать комментарий к своей правке, описывающий проделанные изменения. Также нужно указать источник данных. Активировав флаг Автоматически брать источник из текущих слоёв, можно автоматически заполнить эту строку названием слоёв спутниковых снимков, использованных при оцифровке. Если данные были получены путем обхода, то указываем survey или дописываем через точку с запятой, например, Bing; survey. Нажимаем кнопку Отправить изменения на сервер и наслаждаемся своими изменениями на карте OSM.

Если на одном участке работало несколько картографов, в момент валидации может появиться окно с конфликтом правок, но это уже совсем другая история.

В заключение

Редактирование в JOSM требует некоторого опыта работы с данными OpenStreetMap. Но в большинстве случаев, когда дело доходит до обрисовки большого района, массового исправления тегов, работы с отношениями, JOSM оказывается удобнее редактора iD.

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

Во время последнего локдауна в Великобритании мы с женой играли в GeoGuessr. Эта игра более размеренна, чем те, в которые мы обычно играем, но хорошо подходит для нашей семьи с 11-недельным младенцем, который становится активнее с каждым днём.

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


Нас серьёзно заинтересовали ежедневные соревнования (Daily Challenge) на GeoGuessr. Мы начали заходить на сайт каждый день и пытаться поставить новый рекорд. В формате Daily Challenge на каждый раунд выделяется по три минуты, которые мы тратили или на бешеное кликанье по австралийскому бушу (при этом иногда путая его с Южной Африкой), или на обсуждение того, есть ли в шведском языке буква .

Теперь у меня накопился большой объём знаний типа увижу узнаю. Я могу опознать Гренландию с первого взгляда. Вернулись мои утерянные знания флагов стран, а также появились новые знания о флагах штатов США, о тех странах, где ездят по левой и правой полосам, где используют километры или мили. Я знаю почти все доменные имена стран (их часто можно встретить на рекламных билбордах вдоль дорог) мне ещё долго не забыть .yu.

Вы знали, что чёрно-белые дорожные ограждения распространены в России и Украине? Или что можно разобрать синюю полосу EU на автомобильных номерах, несмотря на размытие Google Street View? Подробнее об этом можно прочитать в этом руководстве из 80 тысяч слов Geoguessr the Top Tips, Tricks and Techniques.

Указывающая вниз полосатая красно-белая стрелка даёт понять, что вы находитесь в Японии, с большой вероятностью, на острове Хоккайдо или, возможно, на острове Хонсю рядом с горами.

Немного глубокого обучения

Однажды я прочитал, что машинное обучение уже умеет делать всё, что и человек, но меньше чем за одну секунду. Распознать лицо, выбрать текст из изображения, повернуть, чтобы не врезаться в другую машину. Это заставило меня задуматься, а размышления привели к статье под названием Geolocation Estimation of Photos using a Hierarchical Model and Scene Classification, написанной Эриком Мюллером-Будаком, Кадером Пусту-Иреном и Ральфом Эвертом. В этой статье геолокализация рассматривается как задача классификации, в которой Земля подразделена на географические ячейки.

Она прогнозирует GPS-координаты фотографий.


Даже по фотографиям, которые сделаны в помещении! (Daily Challenge игры GeoGuessr часто засовывает игрока внутрь музеев).

Недавно авторы статьи выпустили реализацию на PyTorch и указали веса для обученной модели base(M, f*) с внутренней архитектурой ResNet50.

Я предположил, что обученная модель не очень хорошо будет соответствовать тем частям фотосфер, которые я смогу получить от GeoGuessr. В качестве данных обучения авторы использовали подмножество набора данных из 100 миллионов фотографий Yahoo Flickr Creative Commons (YFCC100M). В него вошли примерно пять миллионов изображений Flickr с геометками и неопределённых фотографий, например, снимков внутри помещений, еды и людей, местоположение которых сложно спрогнозировать.

Любопытно было то, что в наборе данных Im2GPS люди определяли местоположение изображения с точностью на уровне страны (в пределах 750 км) в 13,9% случаев, а Individual Scene Networks справлялись с этой задачей в 66,7% случаев!


Итак, возник вопрос: кто лучше в GeoGuessr, моя жена (потрясающий игрок) или машина?

Автоматизируем GeoGuessr с помощью Selenium

Для скрейпинга скриншотов из текущего внутриигрового местоположения я создал программу на Selenium, четыре раза выполняющую следующие действия:

• Сохраняем скриншот canvas
• Делаем шаг вперёд
• Поворачиваем обзор примерно на 90 градусов

Количество повторов этих действий можно настроить через NUMBER_OF_SCREENSHOTS в показанном ниже коде.

'''Given a GeoGuessr map URL (e.g. https://www.geoguessr.com/game/5sXkq4e32OvHU4rf)take a number of screenshots each one step further down the road and rotated ~90 degrees.Usage: "python file_name.py https://www.geoguessr.com/game/5sXkq4e32OvHU4rf"'''from selenium import webdriverimport timeimport sysNUMBER_OF_SCREENSHOTS = 4geo_guessr_map = sys.argv[1]driver = webdriver.Chrome()driver.get(geo_guessr_map)# let JS etc. loadtime.sleep(2)def screenshot_canvas():    '''    Take a screenshot of the streetview canvas.    '''    with open(f'canvas_{int(time.time())}.png', 'xb') as f:        canvas = driver.find_element_by_tag_name('canvas')        f.write(canvas.screenshot_as_png)def rotate_canvas():    '''    Drag and click the <main> elem a few times to rotate us ~90 degrees.    '''    main = driver.find_element_by_tag_name('main')    for _ in range(0, 5):        action = webdriver.common.action_chains.ActionChains(driver)        action.move_to_element(main) \            .click_and_hold(main) \            .move_by_offset(118, 0) \            .release(main) \            .perform()def move_to_next_point():    '''    Click one of the next point arrows, doesn't matter which one    as long as it's the same one for a session of Selenium.    '''    next_point = driver.find_element_by_css_selector('[fill="black"]')    action = webdriver.common.action_chains.ActionChains(driver)    action.click(next_point).perform()for _ in range(0, NUMBER_OF_SCREENSHOTS):    screenshot_canvas()    move_to_next_point()    rotate_canvas()driver.close()

Скриншоты содержат и интерфейс GeoGuessr, но я не стал заниматься его удалением.

Приблизительное определение геолокации

Я перешёл к ветке PyTorch branch, скачал обученную модель и установил зависимости с помощью conda. Мне понравился README репозитория. Раздел requirements был достаточно понятным и на новом Ubuntu 20.04 у меня не возникло никаких проблем.

Для выяснения отношений между человеком и машиной я выбрал в GeoGuessr карту World. Отправив URL своей программе Selenium, я прогнал её для четырёх скриншотов, сделанных в GeoGuessr.

Ниже представлены сокращённые результаты работы машины.

python -m classification.inference --image_dir ../images/                                lat        lngcanvas_1616446493 hierarchy     44.002556  -72.988518canvas_1616446507 hierarchy     46.259434  -119.307884canvas_1616446485 hierarchy     40.592514  -111.940224canvas_1616446500 hierarchy     40.981506  -72.332581

Я показал те же четыре скриншота своей жене. Она предположила, что точка находится в Техасе. На самом деле место находилось в Пенсильвании. Машина сделала для каждого из четырёх скриншотов четыре различные догадки. Все догадки машины находились в США. Две достаточно близко друг к другу и две подальше.


Если взять усреднённое местоположение, то машина в этом раунде побеждает!

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

После написания этого поста я узнал о потрясающей предыдущей работе со сравнением результатов человека и машины на поле боя GeoGuessr. В статье PlaNet Photo Geolocation with Convolutional Neural Networks Тобиас Вейанд, Илья Костиков и Джеймс Филбин пытались определить местоположение фотографии всего по нескольким пикселям.

Чтобы выяснить, насколько PlaNet сравнима с интуицией человека, мы позволили ей соревноваться с десятью много путешествовавшими людьми в игре Geoguessr (www.geoguessr.com).

В сумме люди и PlaNet сыграли в 50 раундов. PlaNet выиграла 28 из 50 раундов с медианной погрешностью локализации в 1131,7 км, в то время как медианная погрешность людей составляла 2320,75 км.

Веб-демо

Авторы статьи Geolocation Estimation of Photos using a Hierarchical Model and Scene Classification создали довольно милый веб-инструмент. Я проверил его на одном из скриншотов Selenium.

Графическая демонстрация, в которой вы сможете посоревноваться против описанной в статье системы с глубоким обучением находится здесь: https://tibhannover.github.io/GeoEstimation/. Также мы создали многофункциональный веб-инструмент, поддерживающий загрузку и анализ пользовательских изображений: https://labs.tib.eu/geoestimation

Обучаемость GeoGuessr

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

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

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

Если в street view посмотреть вниз то можно увидеть часть машины, снимавшей текущую фотосферу. Например, в Кении спереди у машины есть чёрная труба. Основная часть Вьетнама была снята с мотоцикла, и часто можно увидеть шлем водителя. Страны часто снимаются одной машиной с уникальным цветом или антенной.


В других местах в небе есть место, где сшитая фотосфера выглядит разорванной (в основном в Сенегале, Черногории и Албании). В Африке за автомобилем Street View иногда едут машины сопровождения. Есть разные поколения камер, с разным разрешением, типами гало, цветопередачей и размытием в нижней части сферы. В нижнем углу фотосферы есть сообщение об авторстве, обычно это Google и указание года, но иногда бывает и имя фотографа.

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

---

На правах рекламы

Закажите сервер и сразу начинайте работать! Создание VDS любой конфигурации в течение минуты. Эпичненько :)

Всем привет! Меня зовут Константин Измайлов, я руководитель направления Data Science в Delivery Club. Мы работаем над многочисленными интересными и сложными задачами: от формирования классических аналитических отчетов до построения рекомендательных моделей в ленте приложения.

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

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

Статья написана по мотивам выступления с Евгением Макиным на конференции Highload++ Весна 2021. Для тех, кто любит видео, ищите его в конце статьи.

Бизнес-модель работы Delivery Club

Бизнес-модель Delivery Club состоит из двух частей:

• ДДК (доставка Деливери Клаб): мы передаем заказ в ресторан и доставляем его клиенту, то есть ресторану остается только приготовить заказ к определенному времени, когда придет курьер.
• МП (маркетплейс): мы передаем заказ в ресторан, а он своими силами доставляет заказ в пределах своей согласованной зоны.

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

При этом построить зону доставки это полбеды. Ведь география ДДК огромна и охватывает более 170 городов России с тысячами ресторанов в каждом из них, и у каждого ресторана свои индивидуальные особенности, которые необходимо учитывать при построении эффективного сервиса доставки.

Рисуем зону доставки ресторана

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

Как процесс выглядел раньше

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

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



Стоит упомянуть и про SLA (Service Level Agreement соглашение о максимальной длительности отрисовки зоны доставки для одного партнера): онбординг партнера или подготовка его зоны для внедрения в нашу платформу составляли порядка 40 минут для одного заведения. Представьте, что к вам подключилась городская сеть с сотней ресторанов, а если это ещё и жаркий сезон, например, после проведения рекламной акции Вот наглядное доказательство неэффективности ручной отрисовки:

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

Проблемы ручной отрисовки зон:

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

Baseline

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

При этом оставались недостатки:

• артефакты в зонах доставки (стандартный случай с переходом через реку);
• единообразный подход к партнерам (не учитываются индивидуальные KPI партнеров).

Через Delivery Club ежедневно проходят сотни тысяч заказов, и артефакты зон доставки могли дать большую нагрузку на колл-центр.



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

Преимущества технологии H3

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


Источник: eng.uber.com/h3

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

Также стоит отметить, что:

1. Существует хорошая библиотека для работы с H3, которую и выбрала наша команда в качестве основного инструмента. Библиотека поддерживает многие языки программирования (Python, Go и другие), в которых уже реализованы основные функции для работы с гексагонами.
2. Наша реляционная аналитическая база Postgres поддерживает работу с нативными функциями H3.
3. При использовании гексагональной сетки благодаря ряду алгоритмов, работающих с индексами, можно очень быстро получить точную информацию о признаках в соседних ячейках, например, определить вхождение точки в гексагон.
4. Преимуществом гексагонов является возможность хранить информацию о признаках не только в конкретных точках, но и в целых областях.

Алгоритм построения зон доставки

1. Вокруг ресторана строим лучи, вдоль которых создаем набор точек, а затем с помощью маршрутизатора рассчитываем расстояние от ресторана до каждой из точек.
   
   
   
   
   
2. Убираем точки-выбросы. Анализируем все постройки и сразу отбрасываем те, которые нас не интересуют. Например, какие-то мелкие нежилые объекты. Далее с помощью DBSCAN формируем кластеры точек и отбрасываем те, которые для нас не являются важными: например, если кластер находится далеко от ресторана или нам экономически невыгодно доставлять туда.
   
   
   
   
   
3. Далее на основе очищенного набора точек применяем триангуляцию Делоне.
   
   
   
   
   
4. Создаем сетку гексагонов H3.
   
   
   
   
   
5. Определяем, к какому треугольнику принадлежит каждый гексагон H3, и для этого гексагона определяем расстояние от ресторана с помощью аппроксимации по трем вершинам треугольника (кто хочет чуть больше разобраться в геоматематике и понять рациональность использования тех или иных геоформ, вот ссылка на оригинал статьи Uber).
   
   
   
   
   
6. Далее оптимизируем функцию ошибки, которая зависит от нашей задачи, постепенно удаляя с внешней стороны зоны по одному гексагону. При этом мы следим за формой нашей зоны, она должна оставаться цельной.
   
   
   
   
   
   В текущей версии алгоритма мы используем следующие функции ошибок:
   
   
   • минимизацию времени доставки при фиксированном покрытии;
   • максимизацию охвата пользователей при фиксированном времени доставки.
   
   Пример функции ошибки для минимизации времени доставки:
   

   

   
   где функция ошибки минимизации времени доставки с фиксированным покрытием,
   время от ресторана i до клиента,
   количество клиентов в зоне доставки.
   
   
7. Далее строим временной градиент в получившихся зонах (с очищенными выбросами) и с заранее определенными интервалами (например, по 10-минутным отрезкам пешего пути).
   
   
   
   
   

Зона построена, теперь очередь оптимизации бизнес-метрики. В алгоритме используются оптимизации покрытия клиентов с фиксированием времени доставки и оптимизация времени доставки с фиксированным покрытием.

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

Внедрение

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

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

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

Но тут пришел COVID-19

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

Всё это отчасти помогло нам убедить партнёров активнее тестировать и переходить на алгоритм автоматического построения зон доставки. В итоге, за первую неделю после начала пандемии и введения ограничения для ресторанов порядка 50% наших партнеров перешло на новые зоны доставки, которые уже были построены автоматически, а в течение следующих 3 недель уже почти не осталось партнеров с зонами, построенными вручную.

Оценка

После решения всех горящих проблем нам нужно было немного отдышаться и понять, что мы вообще наделали. Для этого воспользовались A/B-тестом, а точнее его вариацией switch-back. Мы сравнивали зоны ресторанов с одинаковыми входными параметрами, оценивали GMV и время доставки, где в качестве контроля у нас были простые автоматически отрисованные зоны в виде окружностей и прямоугольников, либо зоны, отрисованные операторами вручную.



В итоге для разных алгоритмов мы получили прирост GMV в тестовых алгоритмических зонах при одинаковом времени доставки, либо уменьшение времени доставки в случае с алгоритмом с фиксированным клиентским покрытием.

А время, затрачиваемое на построение зон для партнера из примера выше, теперь выглядит более оптимистично:

Ускорение в 670 раз!

Текущая ситуация и планы

Сервис работает в production. Зоны автоматически строятся по кнопке. Появился более гибкий инструмент для работы со стоимостью доставки для клиентов в зависимости от их удаленности от ресторана. 99,9% ресторанов (изредка ещё нужно ручное вмешательство) перешли на алгоритмические зоны, а наш алгоритм поспособствовал переходу бэкенда на H3.

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

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

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

Всем спасибо!