В статье Танцующие горы Ирана по данным спутниковой интерферометрии показан очень необычно выглядящий результат спутниковой интерферометрии. Сегодня мы посмотрим, что же это значит и почему именно этого и следует ожидать. Ранее в статье PyGMTSAR, или спутниковая интерферометрия для всех с примерами Jupyter Python ноутбуков на Google Colab я рассказал про свой пакет для спутниковой интерферометрии на основе радарных снимков Sentinel-1 PyGMTSAR (Python GMTSAR), написанный именно для того, чтобы быстро и удобно получить и проанализировать результаты непосредственно в Python ноутбуке.

По ссылке вы найдете готовый ноутбук на Google Colab, позволяющий прямо в браузере выполнить всю обработку и увидеть результаты и, при желании, тут же поработать с ними: Yamchi DAM Interferograms Persistent Scatterer Interferometry (PSI) Analysis Для Debian Linux я сделал скрипт инициализации облачного инстанса GMTSAR.install.debian10.sh, а на Google Colab ноутбук автоматически установит все необходимые зависимости, просто следуйте подсказкам в ноутбуке.

Введение

В статье PyGMTSAR, или спутниковая интерферометрия для всех с примерами Jupyter Python ноутбуков на Google Colab уже рассказывалось, что собой представляет дифференциальная спутниковая интерферометрия (DInSAR), основанная на обработке пар фазовых изображений радарной космической съемки. Мы будем обсуждать исключительно работу со снимками Sentinel-1, которые за счет очень точно известных орбит космических аппаратов позволяют использовать геометрическое выравнивание для получаемых спутниковых снимков.

Замечу, что получение "живого" примера для задач интерферометрии — задача совсем не тривиальная. В самом деле, для создания показанной ниже анимации необходимо обработать около двух десятков спутниковых снимков размером около 5 ГБ каждая, итого ~100 ГБ "сырых" данных и в разы больше промежуточных результатов при их обработке. Чтобы сделать возможным запуск такой обработки на Google Colab (около 70 ГБ дискового пространства доступно на бесплатном аккаунте), я предварительно скачал все необходимые спутниковые данные и удалил ненужные поляризации (оставив только VV поляризацию из доступных VV, HH и кроссполяризаций) и обрезал изображения — выбрал только sub-swaths 2 (один из трех доступных) и вырезал из него нужные полосы (bursts). Полученные снимки упакованы в один архив и размещены на общедоступном Google Cloud Storage, с которого их и скачивает указанный выше ноутбук. Кроме того, существуют разные подходы к построению интерферометрических пар и мы воспользуемся наименее ресурсоемким, хотя мне и понадобилось предварительно провести вычисления для разных подходов, чтобы оценить необходимые параметры обработки.

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

Более детальное видео доступно на YouTube:

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

PSI подход предполагает использование лишь минимально необходимых пар снимков так, что каждый снимок за исключением мастера используется не более чем в одной паре (PSI), а SBAS основан на построении большого множества возможных пар и автоматизированном их анализе для уменьшения влияния различных помех. Очевидно, второй путь намного более ресурсозатратный, так что мы ограничимся первым. Кстати, в ноутбуке есть подсказки, как можно запустить обработку для обоих методов (я не рекомендую запускать SBAS анализ на бесплатном инстансе Google Colab).

В обоих случаях сначала необходимо с суб-пиксельной точностью выровнять все изображения в единый стэк так, чтобы обеспечить возможность когерентности между пикселами всех изображений. Соответствующее мастер-изображение для выравнивания называется супермастер и располагается в центре поперечной базовой линии (определяемой как расстояние по перпендикуляру между лучом зрения спутника при получении всех используемых изображений), в нашем случае это снимок от 2021-02-13. При ошибках выравнивания интерферограммы получить вовсе не удастся или они будут очень сильно зашумлены и на них окажется возможным разглядеть только элементы рельефа. Мастером или опорным снимком называется первое изображение в каждой паре, и если в технике SBAS каждое изображение должно быть мастером во многих парах, то в PSI мастер один и входит в каждую пару. В нашем случае выбран мастер-снимок с датой 2021-04-26, смотрите SBAS baseline диаграмму ниже:

Метод PSI основан на том, что во многих случаях возможно выбрать один так называемый мастер-снимок и построить интерферометрические пары со всеми оставшимися снимками. При анализе на небольшом интервале времени мастер-снимок может быть первым снимком, тогда мы сразу увидим последовательные изменения на полученных результатах. Под результатами здесь подразумеваются смещения в направлении луча зрения спутника, вычисленные на основании развернутых фаз (unwrapped) построенных интерферограмм. Для увеличении временного интервала уменьшается когерентность между снимками и интерферограммы становятся все более зашумленными. Зачастую, это не проблема — для технологических строений типа дамбы или открытых скалистых участков возможно построить интерферограммы даже по снимкам, сделанным с интервалом в несколько лет. В нашем же случае, когда мы хотим получить картину смещений для большой области, интервал между снимками приходится ограничивать. В таком случае часто используется тоже довольно простой путь с выбором мастер-изображением снимка из середины рассматриваемого временного интервала — в этом случае максимальный временной интервал между любыми двумя снимками примерно равен половине всего интервала наблюдения. Но есть и нюанс — если это изображение окажется зашумлено (атмосферными помехами, к примеру) или просто не очень точно выровнено, качество результата окажется совершенно неудовлетворительным. Это решается смещением на одно или два изображения влево или вправо по временной шкале для выбора наилучшего мастер-изображения так, чтобы получить максимум суммы когерентностей для всех интерферометрических пар. Здесь как раз понадобится провести предварительный SBAS анализ или многократно выполнить PSI для разных мастер-снимков и выбрать наилучший результат.

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

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

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

Обратим внимание, что для мастер-изображения смещение тождественно равно нулю. Тут все еще не просто понять, насколько качественные результаты мы получили для первых и последних интерферограмм. Чтобы оценить точнее, рассмотрим значения когерентности для объекта с ожидаемо высокой когерентностью — например, для плотины Ямчи:

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

Референсная область

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

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

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

Теперь посмотрим на явление приливов твердой Земли (solid Eath tides), вызванных гравитацией Солнца и Луны. Кроме приливов водных эти же силы вызывают смещения и твердых пород, причем амплитуда смещений даже твердых скал может быть порядка одного метра в высоту. В силу того, что интервал съемки не совпадает с интервалом этих приливов, мы наблюдаем множественные эффекты алиасинга с периодами от нескольких недель и до нескольких лет. Использование метода вычитания референсной области исключает это влияние, посколько пространственный масштаб таких приливов велик. И все же, необходимо уделить приливам особое внимание — в случае, когда после вычитания референсной области мы все еще наблюдаем соответствующие эффекты (коррелированные со смещениями референсной области) это означает, что рассматриваемая территория двигается не так, как это ожидается от твердого тела. Иначе говоря, в недрах присутствуют подземные водные резервуары или подземные магматические камеры или и то и другое вместе (можно ориентировочно определить по пространственному масштабу — магматические камеры расположены глубоко и создают эффекты намного большего пространственного масштаба).

Посмотрим графики смещений для дамбы, референсной области и для дамбы за вычетом среднего значения смещений референсной области:

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

Заключение

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

Что здесь можно улучшить? Во-первых, мы рассматриваем LOS смещения, то есть по направлению луча зрения спутника, а не строго вертикальные (хотя большая часть LOS смещения определяется именно вертикальной компонентой в силу геометрии съемки). Чтобы точно выделить именно вертикальные смещения (проекцию вектора), необходимо выполнить рассмотренные выше вычисления для обеих орбит спутников Sentinel-1 (восходящей и нисходящей) и посчитать среднее значение LOS. Впрочем, здесь возникает значимый нюанс — интервал между съемкой с одной орбиты составляет 12 дней и за это время мы видим значительные смещения. Интервал между съемками с двух орбит вдвое меньше, то есть 6 дней, и если за это время происходят значительные смещения (как в нашем случае), потребуется выполнить некоторую интерполяцию, чтобы посчитать значения смещения с интервалом 6 дней для обеих орбит. Зачастую, это выполняется для уже вычисленных значений тренда, поскольку найти качественную аппроксимацию смещений каждой точки поверхности, вызванные как реально происходящими процессами, так и множественным алиасингом — та еще задача.

Кроме того, представляет большой интерес выделение гармонических годовых колебаний для нахождения подземных водных и газонефтяных резервуаров. Зачастую вода представляет собой значительно большую ценность, нежели нефть — например, на пустынных территориях Ирака ищут именно питьевую или пригодную для сельского хозяйства воду, а с нефтью дела там обстоят гораздо лучше, чем с водой. Для этих целей необходимо провести анализ серии данных минимум за два года, выбрав множество мастер-изображений для метода PSI или используя метод SBAS. На бесплатном инстансе Google Colab такую обработку выполнить будет очень затруднительно из-за недостатка как дискового пространства, так и вычислительных ресурсов.

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

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

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