Spectrum

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

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

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

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