После анализа модели Танцующие горы Ирана по данным спутниковой интерферометрии мне захотелось проверить набор гипотез и улучшить качество результатов. Как оказалось, ни один из существующих интерферометрических пакетов не позволяет этого сделать так, как мне нужно. Оценив фронт работ, я решил, что за месяц фулл-тайм работы я смогу написать свою систему спутниковой интерферометрии для радарных снимков Sentinel-1 на основе открытого продукта GMTSAR, реализовав собственные алгоритмы обработки данных и обеспечив удобную работу в среде Jupyter Python. По образованию я радиофизик и мой диплом магистра по моделированию голограмм в оптически нелинейных средах (равно моделированию интерференции) в свое время был признан победителем во всероссийском конкурсе, так что мне удалось уложиться в поставленные сроки и реализовать все запланированное — больше свободного времени на этот проект у меня просто нет. Итак, встречайте PyGMTSAR (Python GMTSAR) — по ссылке вы найдете готовые ноутбуки, которые в один клик можно запустить на Google Colab и прямо в браузере увидеть результаты и, при желании, тут же поработать с ними. Для Debian Linux я сделал скрипт инициализации облачного инстанса GMTSAR.install.debian10.sh, а на Google Colab ноутбуки автоматически установят все необходимые зависимости, что позволяет легко запускать их в "облаках".

Введение

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

Задача разворачивания фазы отнюдь не однозначная и решается разными способами, в том числе, известный продукт SNAPHU — Statistical-Cost, Netowrk-Flow Algorithm for Phase Unwrapping использует методы двумерного растрового роутинга для решения этой задачи. Да, и здесь нужен роутинг, про который я тоже много писал ранее. Будущим методов разворачивания считается трехмерный подход, поскольку таким способом можно получить однозначное решение для всего стэка интерферограмм. Как оказалось, этого же результата можно достичь и с помощью постобработки результатов двумерной развёртки фазы.

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

Small BAseline Subset (SBAS) и Persistent Scatterer Interferometry (PSI)

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

Техника PSI для этого же набора снимков предполагает анализ лишь 4 последовательных интерферограмм. Анализируя все 9 интерферограмм, но выбирая только пикселы с высокой когерентностью для всех пар, мы получаем комбинированный метод SBAS+PSI. Впрочем, следует учитывать, что далеко не все интерферограммы показывают смещения поверхности, поскольку атмосферные помехи могут скрыть все детали поверхности, обратите внимание на субдиагональные полосы на 5ти из 9ти интерферограмм, как показано ниже:

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

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

На самом деле, здесь в стэке снимков лишь один снимок с атмосферными помехами, который и портит 5ть интерферограмм из 9ти. И мы можем его найти методом анализа всех треугольников на SBAS диаграмме. Рассмотрим треугольник из трех снимков A,B,C и трех интерферограмм AB, BC, AC. Очевидно, что сумма изменений за два интервала должна равняться изменению за эти же два интервала, то есть выполняется тождество AB +BC = AC. Если это условие не выполняется, тогда как минимум один из снимков содержит помехи. Составив и решив систему уравнений для всех треугольников на диаграмме SBAS, мы найдем все корректные и зашумленные изображения. Кроме того, можно использовать корреляционный анализ для выбранных триплетов, как показано ниже, в случае зашумленных снимков возникает ложная высокая корреляция одного знака с любым снимком. Смотрите интерферограммы для одного триплета и значения корреляции для всех триплетов, очевидно, здесь зашумленным является снимок 2015-04-27:

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

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

Заключение

В процессе выполнения проекта мной написан программный код для проведения всего необходимого анализа, выполнена проверка на известных примерах и результаты анализа доступны в виде "живых" ноутбуков на Google Colab, смотрите ссылки на странице проекта PyGMTSAR (Python GMTSAR). Кроме того, я реализовал распараллеливание стэкирования снимков и собственные оригинальные алгоритмы, в частности, почти мгновенное матричное преобразование радарных и географических координат, что позволяет на лету выполнять это преобразование в любой момент, как это и сделано в ноутбуках. А еще пример из GMTSAR S1A_Stack_CPGF_T173.ipynb вычисляется втрое быстрее и к нему добавлены средства поиска и исключения зашумленных снимков, см. S1A_Stack_CPGF_T173_TODO.ipynb Обработка рельефа также сделана аккуратно методом интерполяции с контролем непрерывности первой производной, что исключает артефакты, возникающие в GMTSAR (смотрите мои тикеты в багкрекере проекта, там много интересного для специалистов). Многие используемые утилиты из GMTSAR мне пришлось править и процесс приема патчей в основной репозиторий идет, но не быстро, так что с установкой оригинального GMTSAR указанные ноутбуки пока работать не могут.

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

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

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