В дополнение к открытым спутниковым данным, некоторые из которых перечислены в статье Общедоступные данные дистанционного зондирования Земли: как получить и использовать, существует и множество производных продуктов — например, рельеф. Притом можно найти открытый рельеф разного пространственного разрешения, равно как и множество коммерческих, и появляется задача выбрать лучший продукт из доступных.
Сегодня мы рассмотрим открыто доступный рельеф пространственным разрешением 30 м и 1 м и сравним с коммерческим разрешением 1 м. Для сравнения и оценки рельефа разного масштаба используем методы анализа пространственного спектра, неоднократно описанные в моих предыдущих статьях, например, Пространственные спектры и фрактальность рельефа, силы тяжести и снимков В силу фрактальной природы рельефа, его спектр в двойных логарифмических координатах должен совпадать с линией, и мера этого совпадения и есть качество рельефа, а разрешение, начиная с которого компоненты спектра подчиняются закону фрактальности, его реальное разрешение. Как будет показано на высокодетальном открытом рельефе, выбранный метод оценки корректен.
Рельеф USGD NED DEM 1m и ALOS DEM 30m со спутниковой картой Google Satellite
Введение
Занимаясь разными задачами (геологического) моделирования, раз за разом я сталкивался с одними и теми же проблемами в используемом рельефе (как открыто доступном, так и коммерческом). Скажем, гидродинамическое моделирование затопления местности, построение трехмерной геологической модели (решение обратной задачи), спутниковая интерферометрия — все это требует знания детального рельефа местности. Но какой рельеф ни возьми, если детально его рассмотреть, мы увидим огромное количество артефактов. Пусть визуальная оценка это дело субъективное, но и в пространственных спектрах наблюдается удручающая картина.
Рассмотрим ниже некоторые принципиальные особенности цифрового рельефа и далее проверим его качество на примерах от различных производителей.
Заявленное пространственное разрешение рельефа
Пространственное разрешение рельефа, даже построенного по одним и тем же исходным данным, варьируется в широких пределах в зависимости от территории. И, тем более, разное разрешение и точность имеют результаты, полученные разными способами — скажем, методом триангуляции космоснимков и с помощью лидарной съемки. К примеру, для территории США доступно множество продуктов рельефа разрешением от 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 м.
Вертикальная точность получена относительно USGD NED DEM точностью 10 см. Среднеквадратичная ошибка составляет (RMSE) 5.1 м, что полностью соответствует заявленной точности обоих продуктов: 5 м точность рельефа ALOS плюс 0.1 м точность рельефа USGD NED DEM.
Коммерческий рельеф 1 м, построенный методом триангуляции. Реальное разрешение skipped
Образец коммерческого рельефа со спутникового аппарата Pléiades от компании Airbus DS Geo SA получен с официального сайта Elevation1 DSM + Pléiades 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%). Таким образом, как и следует из теории, точный рельеф фрактален на всех масштабах.
Вертикальную точность, очевидно, напрямую проверить мы не можем за отсутствием данных такой детальности. Тем не менее, косвенная проверка — вычисление среднеквадратичного отклонения (RMSE) от рельефа ALOS World 3D — 30m — показывает ожидаемый результат 5.1 м, что полностью соответствует заявленной точности обоих продуктов: 5 м точность рельефа ALOS плюс 0.1 м точность рельефа USGD NED DEM.
Заключение
С помощью качественного рельефа высокого разрешения можно выполнить гидродинамическое моделирование на поверхности, смотрите Гидродинамическое моделирование (CFD) на рельефе с помощью MantaFlow и визуализация результатов в ParaView, и построить детальные геологические модели как показано в статье Построение достоверных геологических моделей. При использовании методов фрактальной математики становится возможным выделение рудных объектов метрового масштаба, см. Пространственные спектры и фрактальность рельефа, силы тяжести и снимков. Также детальный рельеф помогает, при определенных условиях, уточнить спутниковые интерферограммы и получить более детальную картину отражения сейсмических волн от глубинных объектов и более точную модель смещений поверхности, см. Геология XXI века как наука данных о Земле и Вычислительная геология и визуализация.
В следующей статье мы рассмотрим синтез высокодетального рельефа по общедоступным открытым данным, используя статистические методы. В самом деле, нам доступны серии радарных и оптических космических снимков разрешением от 5 м, и улучшение открытого глобального рельефа 30 м до разрешения снимков выглядит очень заманчиво. Вот мы и проверим, как и насколько это возможно сделать.
