Dem

Сразу оговорюсь, что деление на века немного условно. Например, спутниковая интерферометрия используется с конца 1980-х годов, при этом высококачественные данные стали общедоступными только в 2000-х годах. Трехмерные модели тоже отнюдь не новинка, и делали их ну очень давно ведь и плоская Земля на трех китах вполне себе объемная модель. Так в чем же разница геологии века прежнего и настоящего?

Слева фрагмент геологической карты США, справа 3D геологическая модель с интерферограммой на поверхности рельефа по данным радарной спутниковой съемки (на шкале Density Gradient,% является характеристикой неоднородности геологической плотности, а Band Magnitude обозначает разность фаз отраженного сигнала радара для пары разновременных снимков)

Геология: ремесло или наука

Как нам подсказывает википедия:

Геология (от др. -греч. Земля + учение, наука) совокупность наук о строении Земли, её происхождении и развитии, основанных на изучении геологических процессов, вещественного состава, структуры земной коры и литосферы всеми доступными методами с привлечением данных других наук и дисциплин.

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

Век XX

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

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

Век XXI

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

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

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

На следующей картинке показано сечение модели через эпицентр землетрясения правый от центра блок поднялся вверх и левый от центра опустился вниз в результате этого сейсмического события:

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

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

Дисциплина Computational fluid dynamics(CFD) или, на русском языке, Вычислительная гидродинамика изучает поведение различных потоков, в том числе вихревых. Это и моделирование цунами, и лавовых потоков, и выбрасываемых из жерла вулкана камней вместе с лавой и газами и многое другое. Посмотрим, как можно использовать совместно MantaFlow и ParaView, реализовав на встроенном в MantaFlow языке Python необходимые функции конвертации данных. Как обычно, исходный код смотрите в моем GitHub репозитории: MantaFlow-ParaView.

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

Введение

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

Как и в случае с методами спутниковой интерферометрии, см. предыдущую статью Геология XXI века как наука данных о Земле, методы моделирования потоков позволяют узнать многое о происходивших процессах именно узнать, а не предположить. Действительно, моделируя лавовые потоки на палеорельефе (построенным, к примеру, путем решения обратной задачи геофизики), мы можем сравнить полученную модель с реально существующими [застывшими] лавовыми потоками и полями, поскольку они достаточно прочные и могут хорошо сохраняться многие и многие миллионы лет. Построенная модель позволяет изучить залегание различных слоев там, где у нас нет геологической информации, зачастую, один и то же вулкан извергался многократно из разных жерл, при этом потоки лавы разных извержений [и разного состава] могут перекрываться. Кроме того, эта же модель покажет неувязки между используемой моделью палеорельефа и существующими лавовыми проявлениями и позволит внести уточнения. То есть, вместо сложно формализуемой геологической интуиции, мы можем работать с моделью в том числе, изменять параметры и оценивать, насколько разные геологические предположения вообще разумны. Поскольку я сам не геолог, а физик, для меня путь моделирования представляет вполне понятный интерес. В результате, геолог может точнее оценить возможные области залегания полезных ископаемых и их потенциал разумеется, в геологическом исследовании вовсе не ставится задача обойтись без геолога, но, как и везде, ценность результатов согласуется с известным принципом информатики: "Мусор на входе мусор на выходе", поэтому каждая возможность уточнить, заверить и дополнить имеющиеся данные в буквальном смысле на вес золота (или нефти, или воды,...).

MantaFlow

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

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

А это модель заполнения потоком воды заданного рельефа:

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

Добавим рельеф в MantaFlow

Для меня самый удобный вариант это использование ParaView с моим расширением N-Cube ParaView plugin for 3D/4D GIS Data Visualization для построения 3D модели рельефа в ParaView на основе NetCDF или GeoTIFF данных и сохранение нужного участка в формате OBJ для его использования в MantaFlow. Поскольку MantaFlow умеет этот формат загружать и работать с ним в безразмерных координатах, нам потребуется лишь указать нужный размер в безразмерных координатах (скажем, 100% по горизонтальным координатам и 25% по вертикальной чтобы было достаточно места для моделирования столба дыма) и сохранить параметры преобразования для экспорта результатов в физических координатах. Вот скрипт репозитория с реализацией соответствующей функции: mesh2manta.py

Сохраним результаты моделирования в MantaFlow для ParaView

Поскольку мы задаем исходное пространство моделирования в физических координатах (файл OBJ и коэффициенты его масштабирования), у нас есть все необходимое, чтобы и результаты сохранить в физических координатах. По умолчанию, MantaFlow сохраняет безразмерные результаты в формате сжатых массивов Numpy, поэтому мы добавим сохранение в формат с поддержкой физических координат (VTK), см. скрипт репозитория npz2vtk.py. Добавлю, что в скрипте создается массив xarray: N-D labeled arrays and datasets in Python, из которого одной командой можно сохранить данные в формате NetCDF и некоторых других.

Визуализация в ParaView

Как мы уже рассмотрели в предыдущих статьях (с примерами), ParaView поддерживает работу с сериями данных, так что мы можем работать с 4D данными например, в виде 3D анимации. Вот пример анимации вулканического дыма из серии файлов VTK, экспортированных из MantaFlow:

Модели высокой детализации

Увеличение детальности моделей требует и больше ресурсов для их построения. Если модели на сетке 64х64х64 вычисляются за несколько минут на ноутбуке, то при удвоении разрешения по каждой координате время увеличивается в 8 раз (третья степень двойки).

Ниже показана намного более детальная модель на примере турбулентного торнадо с сайта проекта MantaFlow:

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

Заключение

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

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

https://www.youtube.com/playlist?list=PLwr8DnSlIMg0KABru36pg4CvbfkhBofAi

Как-то на Хабре мне попалась довольно любопытная статья Научно-технические мифы, часть 1. Почему летают самолёты?. Статья довольно подробно описывает, какие проблемы возникают при попытке объяснить подъёмную силу крыльев через закон Бернулли или модель подъёмной силы Ньютона (Newtonian lift). И хотя статья предлагает другие объяснения, мне бы всё же хотелось остановиться на модели Ньютона подробнее. Да, модель Ньютона не полна и имеет допущения, но она даёт более точное и интуитивное описание явлений, чем закон Бернулли.

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

Я же решил проверить, что же произойдёт в модели Ньютона если её улучшить. Что если добавить в неё недостающий элемент межатомного взаимодействия? Исходный код и бинарники получившегося симулятора доступны на GitHub.

Computational Fluid Dynamics

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

Если вас пугает внешний вид этих уравнений, но вы хотели бы разобраться, то для вас есть очень хорошее объяснение в 7-м томе лекций Ричарда Фейнмана. Загляните в главы 38 (Упругость), 40 (Течение сухой воды) и 41 (Течение мокрой воды). Если совсем кратко и на пальцах, то система уравнений Навье-Стокса это векторное уравнение второго закона Ньютона. Эта система определяет равнодействующую всех сил (давления, вязкости и гравитации) для всех направлений. Дополнительно может быть задано второе векторное уравнение, обеспечивающее условие несжимаемости, если нужно описать жидкости.

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

Один из многих способов численно решить систему уравнений это применить Метод Конечных Элементов в связке с Адаптивным Сеточным Методом в случае подхода Эйлера. Зачем нужна адаптивность сетки и как её можно реализовать, подробно и доступно рассказали ребята из SpaceX в своём докладе. Эйлеров подход как правило, но не обязательно, применяется для моделирования замкнутых объемов неразрывных сред, т. е. тех, в которых отсутствуют пустые места и включения других сред. Для иных сред, чаще всего реализуют подход Лагранжа, через метод сглаженных частиц (Smoothed-particle hydrodynamics, SPH). Метод активно применяется для моделирования воды с полным набором явлений: брызги, капли, лужи, смачивание поверхностей и т. д. Можно даже сымитировать пену или пузыри, если включить в модель частицы воздуха. Реконструкцию поверхности, а точнее изоповерхности, можно произвести любым интересующим вас способом (screen-space meshes, dual contouring, marching tetrahedra, metaballs). Если вы знаете другие интересные подходы, добро пожаловать в комментарии.

Discrete Element Method

Моей задачей было найти такой подход, который позволит обсчитывать миллионы частиц и наблюдать за эволюцией системы практически в реальном времени.

Метод дискретного элемента (Discrete Element Method, DEM) показался мне очень привлекательным, потому что он решает практически схожую задачу. Однако он в первую очередь предназначен для сыпучих и гранулированных сред, где у каждой частицы в трёхмерном пространстве есть позиция, ориентация, произвольная форма и масса.

Чтобы упростить вычисления, я решил оставить у частиц только две степени свободы (координаты X и Y), одинаковую массу и радиус. При построении своей модели в угоду производительности я хотел отбросить параметры и факторы, которые порождают эффекты второго порядка. Однако при моделировании сложных систем, они могут быть очень существенны. Один из показательных примеров это использование NASA модели идеального газа вместо реального при проектировании космических челноков. В результате, во время миссии STS-1 проявились различные аномалии при входе в атмосферу. Подробнее в разделе Mission Anomalies.

Тем не менее у DEM есть одна важная особенность это обнаружение столкновений постфактум (Discrete Collision Detection). Разрешение столкновений происходит простым силовым воздействием по закону Гука.

В противоположность этому подходу существует априорный метод Continuous Collision Detection (CCD), который рассчитывает, когда столкновение произойдёт в будущем. Зная точное время контакта, можно скорректировать временной шаг, и избежать неприятных физических артефактов. Метод активно применяется в современных играх. Для игр CCD очень важен чтобы объекты не туннелировали друг через друга, не проваливались друг в друга и не застревали в текстурах. Метод поддерживается современными движками, в Unity и в Unreal точно.

Подробный доклад о методе Continuous Collision Detection

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

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

Искомый эффект пушечного ядра

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

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

Ударная волна создаёт щит (изображение National Geographic)

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

Подробнее по ссылке с таймкодом https://youtu.be/cx8XbaQNnxw?t=2206

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

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

Температурная карта давления (скалярная сумма модулей сил)

Та же карта, только при большем масштабе

Температурная карта ускорений частиц

Архитектура симулятора

Давайте теперь рассмотрим, как устроен симулятор. Приложение состоит из двух независимых модулей:

1. Рендерер текущего состояния.
2. Модуль симуляции на CPU или CUDA.

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

Фазовое пространство

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

Такое представление очень удобно, потому что практически любую симуляцию можно выразить в таком виде. Движение точки состояния как правило выражается через обыкновенное дифференциальное уравнение первого порядка (Ordinary Differential Equation, ODE). Это уравнение имеет вид $inline$dx/dt = f(x, t)$inline$, где $inline$x$inline$ это позиция точки в фазовом пространстве, а $inline$f$inline$ чёрный ящик, способный определить скорость изменения состояния. Зная $inline$x_0$inline$ и $inline$dx/dt$inline$, можно посчитать следующее значение $inline$x_1 = x_0 + \frac{dx}{dt}dt = x_0 + dx$inline$.

Вне зависимости от того, как работает f, симуляция это всего лишь процесс численного интегрирования.

Подробнее по теме фазового пространства можно ознакомиться в разделах 'Differential Equation Basics' and 'Particle Dynamics' курса https://www.cs.cmu.edu/~baraff/sigcourse/

На канале 3Blue1Brown также доступны отличные материалы:
https://www.youtube.com/playlist?list=PLZHQObOWTQDNPOjrT6KVlfJuKtYTftqH6

Интегратор

После различных экспериментов я решил остановиться на самом грубом, но в тоже время самом простом методе Эйлера (Forward Euler). Я пробовал использовать метод Рунге-Кутты 4-го порядка (RK4), в том числе и с адаптивным шагом, но для конкретно этого сценария больше подошёл метод Эйлера. Преимущество RK4 в том, что он позволяет делать огромные временные шаги ценой четырёхкратного увеличения вычислений, что в некоторых сценариях оправданно. В моём же случае оказалось, что я привязан к малым временным шагам, из-за необходимости избегать туннелирования частиц друг через друга. Кстати, как работают интеграторы с адаптивным временным шагом опираясь на ошибку, можно почитать в 'Differential Equation Basics' lecture notes, section 3, 'Adaptive Stepsizes'.

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

float CSimulationCpu::ComputeMinDeltaTime(float requestedDt) const{    auto rad = m_state.particleRad;    auto velBegin = m_curOdeState.cbegin() + m_state.particles;    auto velEnd = m_curOdeState.cend();    return std::transform_reduce(std::execution::par_unseq, velBegin, velEnd, requestedDt, [](const auto t1, const auto t2)    {        return std::min(t1, t2);    }, [&](const auto& v)    {        auto vel = glm::length(v);        auto radDt = rad / vel;        return radDt;    });}float CSimulationCpu::Update(float dt){    dt = ComputeMinDeltaTime(dt);    m_odeSolver->NextState(m_curOdeState, dt, m_nextOdeState);    ColorParticles(dt);    m_nextOdeState.swap(m_curOdeState);    return dt;}

Вычисление производной состояния

Теперь перейдём к сердцу симулятора определению той самой функции f, упомянутой в параграфе Фазовое пространство. Ниже приведён высокоуровневый код солверов производной для CPU и CUDA версий. Стоит отметить, что CPU версия исторически появилась раньше, так как на ней было проще отладить математику. В CUDA версии появились некоторые улучшения и оптимизации, но суть осталась та же. Отличие состоит в переупорядочивании частиц. Подробнее в разделе Реордеринг частиц.

//CPU-версия void CDerivativeSolver::Derive(const OdeState_t& curState, OdeState_t& outDerivative) {    ResetForces();    BuildParticlesTree(curState);    ResolveParticleParticleCollisions(curState);    ResolveParticleWingCollisions(curState);    ParticleToWall(curState);     ApplyGravity();    BuildDerivative(curState, outDerivative);} //CUDA-версия void CDerivativeSolver::Derive(const OdeState_t& curState, OdeState_t& outDerivative) {     BuildParticlesTree(curState);    ReorderParticles(curState);    ResetParticlesState();    ResolveParticleParticleCollisions();    ResolveParticleWingCollisions();    ParticleToWall();    ApplyGravity();    BuildDerivative(curState, outDerivative);}

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

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

Один из возможных вариантов это использование Uniform Grid, то есть однородной сетки из ячеек на подобии шахматной доски. Одна из реализаций для GPU описана в статье Chapter 32. Broad-Phase Collision Detection with CUDA.

Каждая ячейка пространства содержит в себе список объектов (изображение Tero Karras, NVIDIA Corporation)

В этом случае, поиск столкновений в среднем будет занимать порядка $inline$O(1)$inline$. Каждой частице нужно обойти списки в 9 (3x3) или 27 (3x3x3) ячейках для 2D или 3D случая соответственно. Ещё один приятный плюс структуры это относительная простота распараллеливания её построения. Память под списки можно выделить либо заранее в виде массива, и вычислять выходной индекс через атомарный инкремент, либо строить классический RCU lock-free односвязный список. Nvidia в своих видеокартах уже давно добавила поддержку кучи, поэтому можно вызвать malloc()/free() прямо в device коде, выделяя и освобождая элементы списков.

CppCon 2017: Fedor Pikus Read, Copy, Update, then what? RCU for non-kernel programmers

Однако, у этой структуры есть следующий ряд фундаментальных ограничений:

1. Множество значений координат ограничено размером самой сетки.
2. Близкие ячейки в евклидовом пространстве как правило расположены далеко в адресном пространстве RAM/VRAM, не разделяя единую кэш-линию, что создаёт дополнительную нагрузку на шину памяти.
3. При низкой плотности объектов или малом их количестве структура данных начинает потреблять больше памяти, чем сами данные.
4. Возможно появление чрезмерно длинных списков при большой плотности объектов.
5. В связи с аппаратными особенностями планирования потоков на GPU, некоторые lock-free структуры не способны работать корректно (https://youtu.be/86seb-iZCnI?t=2311, ссылка с таймкодом).

Другой вариант, который я решил использовать в этой симуляции это BVH-дерево на основе Z-кривой. Я наткнулся на эту довольно любопытную структуру данных, когда искал альтернативы однородной сетке.

Первая важная особенность этой структуры данных в её основе лежит фрактальная Z-кривая, она же Кривая Мортона.

Фрактальная Z-кривая (изображение Wikipedia)

Принцип вычисления индекса на кривой чередование битов координат
(изображение Wikipedia)

Задача этой кривой, как и любой другой space-filling curve, состоит в том, чтобы упаковать пространства высших размерностей в одномерное пространство. Если присвоить каждому объекту в 2D/3D пространстве индекс на любой такой кривой, а затем отсортировать все объекты по этому индексу, то мы увидим, что объекты, расположенные близко в геометрическом пространстве, как правило будут лежать близко и в одномерном пространстве. Это свойство позволяет существенно снизить нагрузку на шину памяти. Кстати, если вам нужно обрабатывать изображения, выполняя различные свёрточные операции и применяя фильтры, возможно, вам стоит хранить пиксели в виде одной из такой кривых, а не в виде матрицы.

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

Детально алгоритм описан в статье Maximizing Parallelism in the Construction of BVHs, Octrees, and k-d Trees.
Краткое изложение:

1. Обход дерева: https://developer.nvidia.com/blog/thinking-parallel-part-ii-tree-traversal-gpu/
2. Построение дерева: https://developer.nvidia.com/blog/thinking-parallel-part-iii-tree-construction-gpu/

В сухом остатке, алгоритм следующий. Для каждого из N объектов запускается отдельный поток, который вычисляет bounding box и, на основе его центра, вычисляет код Мортона (индекс на Z-кривой). После этого этапа боксы сортируются в порядке возрастания кода.

Формирование кодов Мортона (изображение Tero Karras, NVIDIA Corporation)

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

Последовательный алгоритм построения префиксного дерева (изображение Tero Karras, NVIDIA Corporation)

Параллельный алгоритм построения префиксного дерева (изображение Tero Karras, NVIDIA Corporation)

После того, как связи между узлами выстроены, начинается этап формирования BVH-структуры.

BVH-структура (изображение Tero Karras, NVIDIA Corporation)

N потоков стартуют с листьев и, поднимаясь к корню, обновляют боксы промежуточных узлов. Так как не определено, какой из детей придёт к родителю первым, то в промежуточных узлах хранится специальный флаг, изначально установленный в ноль. Оба ребёнка с помощью атомарной функции atomicExch() устанавливают флаг в 1. Функция возвращает старое значение, которое было до модификации. Если ребёнку функция вернула 0, то значит он первый. Это также означает, что текущему потоку нельзя модифицировать бокс родителя, потому что бокс его сиблинга может быть ещё не готов. На этом этапе поток завершает своё исполнение. Если же ребёнку функция вернула 1, то можно смело модифицировать родительский бокс, объединяя боксы обоих сиблингов, и снова повторить процесс.

После этого этапа дерево готово к осуществлению запросов.

Реакция на столкновения

В симуляции существует два типа столкновений частица-частица и частица-сегмент профиля.

Реакция частица-частица использует факт того, все объекты уже сохранены в дереве, поэтому существует частная процедура рефлексивного обхода, когда листья ищут столкновения друг с другом. Эта оптимизация была предложена Tero Karras. Особенность процедуры в том, что она распознаёт столкновения A-B и B-A как одно и то же столкновение, поэтому оно детектируется только один раз. Для этого при построении дерева вводится дополнительная информация. В промежуточных узлах хранится индекс самого правого листа (rightmost leaf), до которого можно добраться. Например, на рисунке выше rightmost(N2) = 4, а rightmost(N3) = 8. Когда поток, связанный с листом, скажем, O6, будет опускаться от корня, он обратится к промежуточному узлу N2. Благодаря переменной rightmost, он увидит, что лист O6 недостижим из поддерева N2. В этом случае поток O6 должен проигнорировать всё поддерево N2. Однако, потоки, связанные с листьями из поддерева N2, будут проверять поддерево N3. В конечном итоге, если столкновение с O6 и существует, то об этом сообщит только один поток, и он будет из поддерева N2.

Для частной процедуры рефлексивного обхода прототип функции выглядит следующим образом:

template<typename TDeviceCollisionResponseSolver, size_t kTreeStackSize>void CMortonTree::TraverseReflexive(const TDeviceCollisionResponseSolver& solver);

Для случая частица-сегмент профиля, используется универсальная версия:

template<typename TDeviceCollisionResponseSolver, size_t kTreeStackSize>void CMortonTree::Traverse(const thrust::device_vector<SBoundingBox>& objects, const TDeviceCollisionResponseSolver& solver);

Здесь TDeviceCollisionResponseSolver это объект, который должен реализовать следующий интерфейс:

struct Solver{    struct SDeviceSideSolver    {        ...         __device__ SDeviceSideSolver(...);        __device__ void OnPreTraversal(TIndex curId);        __device__ void OnCollisionDetected(TIndex leafId);        __device__ void OnPostTraversal();    };    Solver(...);    __device__ SDeviceSideSolver Create();}; 

Для каждого тестируемого на столкновение объекта, или листа в случае рефлексивного подхода, создаётся отдельный поток. Каждый поток создаёт свой солвер через фабричную функцию Create(). Далее вызывается метод OnPreTraversal, куда передаётся индекс тестируемого объекта. Если бокс текущего тестируемого объекта перекрыл бокс какого-то листа, вызывается функция OnCollisionDetected с индексом листа. Эта функция отвечает за расчёт физики. После обхода дерева вызывается OnPostTraversal.

Такой формат разрешения коллизий появился неслучайно. С самого начала я реализовал его по-другому. Я разделил обход дерева и вычисление физики на две различные стадии, как это сделал Tero Karras. Однако я столкнулся с проблемой построения списков найденных столкновений. Я попробовал сохранять информацию о коллизиях в виде матрицы NxO, где N количество тестируемых объектов, O максимальный размер списка. Но я отказался от этой идеи, потому что при определенных сценариях быстро заканчивалось место в списках. А это в свою очередь создавало различные физические артефакты. К тому же я обратил внимание, что профилировщик сигнализировал о неэффективной работе с памятью (coalesced memory access). Поэтому я решил попробовать подход без списков, который был описан выше. К моему удивлению, способ оказался немного быстрее и без артефактов.

struct SParticleParticleCollisionSolver{    struct SDeviceSideSolver    {        CDerivativeSolver::SIntermediateSimState& simState;        TIndex curIdx;        float2 pos1;        float2 vel1;        float2 totalForce;        float totalPressure;        __device__ SDeviceSideSolver(CDerivativeSolver::SIntermediateSimState& state) : simState(state)        {        }        __device__ void OnPreTraversal(TIndex curLeafIdx)        {            curIdx = curLeafIdx;            pos1 = simState.pos[curLeafIdx];            vel1 = simState.vel[curLeafIdx];            totalForce = make_float2(0.0f);            totalPressure = 0.0f;        }        __device__ void OnCollisionDetected(TIndex anotherLeafIdx)        {            const auto pos2 = simState.pos[anotherLeafIdx];            const auto deltaPos = pos2 - pos1;            const auto distanceSq = dot(deltaPos, deltaPos);            if (distanceSq > simState.diameterSq || distanceSq < 1e-8f)                return;            const auto vel2 = simState.vel[anotherLeafIdx];            auto dist = sqrtf(distanceSq);            auto dir = deltaPos / dist;            auto springLen = simState.diameter - dist;            auto force = SpringDamper(dir, vel1, vel2, springLen);            auto pressure = length(force);            totalForce += force;            totalPressure += pressure;            atomicAdd(&simState.force[anotherLeafIdx].x, -force.x);            atomicAdd(&simState.force[anotherLeafIdx].y, -force.y);            atomicAdd(&simState.pressure[anotherLeafIdx], pressure);        }        __device__ void OnPostTraversal()        {            atomicAdd(&simState.force[curIdx].x, totalForce.x);            atomicAdd(&simState.force[curIdx].y, totalForce.y);            atomicAdd(&simState.pressure[curIdx], totalPressure);        }    };    CDerivativeSolver::SIntermediateSimState simState;    SParticleParticleCollisionSolver(const CDerivativeSolver::SIntermediateSimState& state) : simState(state)    {    }    __device__ SDeviceSideSolver Create()    {        return SDeviceSideSolver(simState);    }};void CDerivativeSolver::ResolveParticleParticleCollisions(){    m_particlesTree.TraverseReflexive<SParticleParticleCollisionSolver, 24>(SParticleParticleCollisionSolver(m_particles.GetSimState()));    CudaCheckError();}

Во время отладки я обратил внимание, что при высокой плотности частиц, функция OnCollistionDetected как правило вызывается для одних и тех же аргументов среди потоков одного варпа. Типовой сценарий был следующий: если в какой-то области пространства есть частицы A, B, C и D, которые в указанном порядке расположены на Z кривой, то приблизительно происходило вот что:

Как видно из таблицы, на шаге 1 и 2 потоки #1 и #2 выполняли атомарные обращения atomicAdd с одними и тем же частицам C и D в процессе работы функции OnCollistionDetected. Это создаёт дополнительную нагрузку на atomic транзакции.

Начиная с архитектуры Volta, Nvidia добавила в чипы поддержку новых warp-vote инструкций. С помощью инструкции match_any поток может опросить весь warp, получив битовую маску потоков, у которых значение запрашиваемой переменной имеет такое же значение.

Результат работы match_any и match_all для двух кооперативных групп

Коммуникация внутри варпа тоже стала удобнее, потому что появились новые warp shuffle функции с поддержкой маски потоков.

Warp-wide редукция с помощью старых функций без маски

Благодаря этим функциям, потоки перед обращением в глобальную память могут сгруппироваться по признаку общего выходного адреса. Далее эта группа должна выполнить суммирование на уровне регистров SM и уже после этого только один поток обращается в глобальную память. К сожалению, на моём домашнем Pascal (1080 Ti) таких инструкций нет, поэтому я решил попробовать их проэмулировать. Увы, никакого прироста, как и замедления это не дало. Профилировка показала, что хоть нагрузка на atomic транзакции и упала в несколько раз, существенно возросла нагрузка на Arithmetic Workload и увеличилось количество регистров на поток. Заняться разработкой на чипах с Volta или Turing пока не представилось возможным. Хотя, мне всё же удалось протестировать симуляцию на RTX 2060 и найти баг связанной с atomic операцией. Об этом в разделе Барьер памяти.

Другие оптимизации и дополнения

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

SoA

Structure of Arrays одна из техник, которая позволяет ускорить доступ к памяти в определённых ситуациях.

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

typedef uint32_t TIndex; struct STreeNodeSoA {    const TIndex leafs;    int* __restrict__ atomicVisits;     TIndex* __restrict__ parents;     TIndex* __restrict__ lefts;     TIndex* __restrict__ rights;     TIndex* __restrict__ rightmosts;     SBoundingBox* __restrict__ boxes;     const uint32_t* __restrict__ sortedMortonCodes; };

Тоже самое касается и внутреннего состояния солвера производной состояния:

struct SIntermediateSimState {     const TIndex particles;     const float particleRad;     const float diameter;     const float diameterSq;     float2* __restrict__ pos;     float2* __restrict__ vel;     float2* __restrict__ force;     float* __restrict__ pressure; }; 

В тоже время, массив bounding boxов нет смысла хранить в SoA стиле, потому что во всех сценариях необходимы обе точки. Поэтому боксы хранятся в виде Array of Structures (AoS):

struct SBoundingBox {     float2 min;     float2 max; }; struct SBoundingBoxesAoS {     const size_t  count;     SBoundingBox* __restrict__ boxes; }; 

Реордеринг частиц

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

Такой паттерн доступа очень медленно работает на GPU. Для восстановления coalesced memory access, позиции и скорости частиц тоже упорядочиваются вдоль кривой. А после выполнения всех расчётов, силы и давления как выходные величины возвращаются к исходному порядку. Идея не нова и была позаимствована из уже упомянутого доклада SpaceX: https://youtu.be/vYA0f6R5KAI?t=1939 (ссылка с таймкодом).

Восстановление объединённого доступа к памяти (изображение SpaceX)

Такая оптимизация даёт 50% прироста производительности: с 8 FPS до 12 FPS для двух миллионов частиц.

Стек в Shared Memory

Оригинальная статья приводит пример реализации, где стек для обхода дерева реализуется в виде локального массива в скоупе функции. Проблема этого подхода в том, что задействуется локальная память потока область в глобальной памяти. А значит SM начинает выполнять долгие транзакции на чтение и запись, которые ко всему прочему могут оказаться ещё uncoalesced. Суть данной оптимизации, чтобы использовать сверхбыструю Shared Memory на кристалле самого Streaming Multiprocessorа.

__device__ void traverseIterative( CollisionList& list,                                   BVH& bvh,                                    AABB& queryAABB,                                    int queryObjectIdx){    // Allocate traversal stack from thread-local memory,    // and push NULL to indicate that there are no postponed nodes.    NodePtr stack[64]; //<---------------------------- Проблемное место    NodePtr* stackPtr = stack;    *stackPtr++ = NULL; // push    // Traverse nodes starting from the root.    NodePtr node = bvh.getRoot();    do    {        // Check each child node for overlap.        NodePtr childL = bvh.getLeftChild(node);        NodePtr childR = bvh.getRightChild(node);        bool overlapL = ( checkOverlap(queryAABB,                                        bvh.getAABB(childL)) );        bool overlapR = ( checkOverlap(queryAABB,                                        bvh.getAABB(childR)) );        // Query overlaps a leaf node => report collision.        if (overlapL && bvh.isLeaf(childL))            list.add(queryObjectIdx, bvh.getObjectIdx(childL));        if (overlapR && bvh.isLeaf(childR))            list.add(queryObjectIdx, bvh.getObjectIdx(childR));        // Query overlaps an internal node => traverse.        bool traverseL = (overlapL && !bvh.isLeaf(childL));        bool traverseR = (overlapR && !bvh.isLeaf(childR));        if (!traverseL && !traverseR)            node = *--stackPtr; // pop        else        {            node = (traverseL) ? childL : childR;            if (traverseL && traverseR)                *stackPtr++ = childR; // push        }    }    while (node != NULL);}

template<typename TDeviceCollisionResponseSolver, size_t kTreeStackSize, size_t kWarpSize = 32>__global__ void TraverseMortonTreeKernel(const CMortonTree::STreeNodeSoA treeInfo, const SBoundingBoxesAoS externalObjects, TDeviceCollisionResponseSolver solver){    const auto threadId = blockIdx.x * blockDim.x + threadIdx.x;    if (threadId >= externalObjects.count)        return;    const auto objectBox = externalObjects.boxes[threadId];    const auto internalCount = treeInfo.leafs - 1;    __shared__ TIndex stack[kTreeStackSize][kWarpSize]; //Тот самый стек    unsigned top = 0;    stack[top][threadIdx.x] = 0;    auto deviceSideSolver = solver.Create();    deviceSideSolver.OnPreTraversal(threadId);    while (top < kTreeStackSize) //top == -1 also covered    {        auto cur = stack[top--][threadIdx.x];        if (!treeInfo.boxes[cur].Overlaps(objectBox))            continue;        if (cur < internalCount)        {            stack[++top][threadIdx.x] = treeInfo.lefts[cur];            if (top + 1 < kTreeStackSize)                stack[++top][threadIdx.x] = treeInfo.rights[cur];            else                printf("stack size exceeded\n");            continue;        }        deviceSideSolver.OnCollisionDetected(cur - internalCount);    }    deviceSideSolver.OnPostTraversal();}

Использование Shared Memory позволяет достичь прироста на 43%: с 14 FPS до 20 FPS. Подробнее о доступных типах памяти можно почитать в официальной документации:

https://docs.nvidia.com/cuda/cuda-c-programming-guide/index.html#device-memory-accesses

Барьер памяти

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

Артефакт на 20-й RTX серии. Позиции и размер артефактов каждый шаг менялись.

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

Доклад об атомиках и барьерах памяти.

Половина доклада посвящена идеи барьеров памяти и почему они важны при работе с atomic-операциями и lock-free структурами. Дело в том, что процессоры имеют тенденцию переупорядочивать выполнение инструкций (Out-of-order execution), но при этом отслеживая и сохраняя зависимости между ними, чтобы гарантировать корректность. В случае с lock-free структурами для процессоров зависимость не очевидна. Поэтому, нужны барьеры памяти, которые подсказывают процессору, что инструкции не могут быть переупорядочены через барьер. Каждая платформа реализует барьеры по-своему, но, к счастью, разработчики стандарта C++ смогли построить наиболее общую модель. Подробное описание каждой семантики барьеров доступно в документации по std::memory_order.

__device__ void CMortonTree::STreeNodeSoA::BottomToTopInitialization(size_t leafId){    auto cur = leafs - 1 + leafId;    auto curBox = boxes[cur];    while (cur != 0)    {        auto parent = parents[cur];        //1. Опасная atomic операция        auto visited = atomicExch(&atomicVisits[parent], 1);        if (!visited)            return;        TIndex siblingIndex;        SBoundingBox siblingBox;        TIndex rightmostIndex;        TIndex rightmostChild;        auto leftParentChild = lefts[parent];        if (leftParentChild == cur)        {            auto rightParentChild = rights[parent];            siblingIndex = rightParentChild;            rightmostIndex = rightParentChild;        }        else        {            siblingIndex = leftParentChild;            rightmostIndex = cur;        }        siblingBox = boxes[siblingIndex];        rightmostChild = rightmosts[rightmostIndex];        SBoundingBox parentBox = SBoundingBox::ExpandBox(curBox, siblingBox);        boxes[parent] = parentBox;        rightmosts[parent] = rightmostChild;        cur = parent;        curBox = parentBox;        //2. Спасительный барьер памяти.         //Следующая итерация гарантированно увидит результаты всех записей         __threadfence();    }}

Моя ошибка была в том, что я не использовал никаких барьеров памяти в коде, который строит BVH дерево, но при этом активно использует атомарный флаг. Интересно, что оригинальная статья также не использует никаких барьеров. Скорее всего, помимо новой SIMT модели (разделы Volta SIMT Model и Starvation-Free Algorithms) Nvidia добавила в новые архитектуры начиная с Volta более агрессивную реализацию уже упомянутой Out-of-order execution. Как следствие, операции, которые должны были выполняться до atomicExch(), т.е. ещё на предыдущей итерации цикла, на Turing исполняются уже после. В результате такого агрессивного реордеринга инструкций, второй ребёнок, приходя к родителю думает, что его сиблинг уже вычислил и сохранил бокс в общую память, а на самом деле нет. В результате получается data race.

thrust::device_vector

Я слишком поздно заметил, что thurst::device_vector работает в синхронном режиме. Этот контейнер в своём конструкторе и в методе resize() выполняет полную синхронизацию с GPU через cudaDeviceSynchronize(). Видимо, разработчики руководствовались следующими рассуждениями. Раз вектор на видеокарте, то и конструкторы элементов нужно тоже вызывать на видеокарте. Но так как конструкторы могут кидать исключения, нужно дождаться их исполнения, чтобы словить все исключения. Единственный доступный способ для них полная синхронизация. Ещё одна из обнаруженных проблем редукция (reduce, sum, min, max) тоже синхронная, так как возвращает результат на хост. Библиотека Cuda UnBound (CUB) в этом плане куда продуманнее. Кстати, недавно она вошла в состав thrust как бэкенд, хотя раньше её нужно было скачивать отдельно.

Результаты профилировки

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

Картинка кликабельна, можно посмотреть в высоком разрешении

Заключение

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

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

Если вы хотели бы начать изучать CUDA, но не знаете, с чего начать, на Youtube есть отличный курс от Udacity Intro to Parallel Programming. Рекомендую к ознакомлению.
https://www.youtube.com/playlist?list=PLAwxTw4SYaPnFKojVQrmyOGFCqHTxfdv2

На последок, ещё несколько видео симуляций:

CPU-версия, 8 потоков, 131'072 частиц

CUDA-версия, 4.2М частиц, 30 минут симуляции

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

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

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

Существует множество способов определения фрактальной размерности, и общего у них только наличие характерного пространственного масштаба. Можно производить вычисления фрактальности в пространственной области, а можно и в частотной (компоненты пространственного спектра), дифференцировать и интегрировать Секрет в том, что всё многообразие способов дает схожие результаты и выбор конкретного способа довольно произволен, важно лишь понимать, как результаты выбранного метода вычислений сопоставить с результатами других методов (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)