2. sci
  3. Построение карты вертикальных смещений по паре радарных снимков КА Sentinel-1

Построение карты вертикальных смещений по паре радарных снимков КА Sentinel-1

Date Категория sci Теги InSAR / SNAP

Источник данных

Выбор интересующего региона из KML-файла

  1. Выбираем координаты из элемента <coordinates>.
  2. Делаем замены, чтобы привести координаты полигона к виду:
    • Copernicus: footprint:"Intersects(POLYGON((-4.53 29.85, 26.75 29.85, 26.75 46.80,-4.53 46.80,-4.53 29.85)))". Координаты первой и последней точек должны совпадать.
    • Vertex: -4.53,29.85,26.75,29.85,26.75,46.80,-4.53,46.80,-4.53,29.85. Координаты первой и последней точек должны совпадать.
  3. Копируем полученную строку в поле поискового запроса сервиса.
  4. Указываем параметры поиска (начальную и конечную даты, спутник, тип продукта, тип инструмента) — и получаем результаты.

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

  1. Пространственная базовая линия (Bperp) не должна быть выше определенного критического значения (для Sentinel-1 IW — около 200 м). Как правило, этому требованию удовлетворяют снимки, относящиеся к одной относительной орбите (Relative Orbit). Каждый из спутников SENTINEL-1 "наматывает" 175 таких орбит за 12-дневный цикл облета планеты.
    • Узнать номер относительной орбиты в результатах поиска Copernicus: View Product Details > Product > Relative Orbit
  2. Снимки должны относится к одному виду движения спутника: восходящему или нисходящему (ascending / descending).
  3. Снимки не должны отстоять друг от друга на слишком большой промежуток времени. Иначе происходит временнАя декорреляция снимков.

В Vertex существует ASF Interferometry Baseline Tool, позволяющий проверить указанные величины до скачивания снимков.

Проверить длину пространственной (Bperp) и временнОй (Btemp) базовых линий в SNAP: Radar > Interferometric > InSAR Stack Overview, добавить снимки, нажать Overview. Получаем Bperp в метрах и Btemp в сутках.

Итог: нужно выбирать снимки, относящиеся к одной относительной орбите, и желательно отстоящие друг от друга не более чем на 12 суток.

Этапы построения карты смещений на основе снимков IW

  1. Корегистрация.
  2. Генерация интерферограммы.
  3. TOPS deburst.
  4. Выделение интересующей области. <Необязательно>
  5. Topo phase removal.
  6. Multilooking. <Необязательно>
  7. Goldstein filtering.
  8. Развертка фазы интерферограммы.
  9. Преобразование фазы в смещение.
  10. Геокодирование.
  11. Экспорт результатов.

1. Корегистрация

  1. Открыть файлы данных: File > Open Product
  2. Выполнить корегистрацию: Radar > Coregistration > S1 TOPS Coregistration > S1 TOPS Coregistration
    • Выбрать нужные снимки (Source Product) — вкладки Read, Read(2)
    • Выбрать subswath и burst'ы. Используйте для этого встроенную мини-карту. subswath'ы и burst'ы для обоих снимков должны быть одинаковы.

Выбор фрагмента из TOPS-продуктов

Для корегистрации и формирования интерферограммы TOPS-снимки должны обрабатываться "burst by burst". Команда TOPS Split позволяет выбрать: 1) swath, 2) поляризацию и 3) burst'ы. Используйте выбор нужного swath'а и burst'ов вместо задания фрагмента снимка.

Выбрать нужные burst'ы можно c помощью карты в World View (удобно предварительно посмотреть интересующий фрагмент в Google Earth). Запомнив расположение участка, ищем его на мини-карте S1 TOPS Coregistration, вкладка TOPSAR Split. Выбираем burst'ы, перетаскивая границы горизонтального слайдера.

Рекомендуется выбрать поляризацию VV (вкладка TOPSAR Split), см. ошибку Snaphu Import.

По окончании S1 TOPS Coregistration снова смотрим World View и убеждаемся, что выбрали burst'ы, которые содержат интересующий фрагмент снимка.

Шаги 1-3 (Корегистрация > Генерация интерферограммы > TOPS deburst) удобно запускать пакетом: Graph Builder > Load > TOPSAR Coreg Interferogram.

2. Генерация интерферограммы

Radar > Interferometric > Products > Interferogram Formation

3. TOPS deburst

Radar > Sentinel-1 TOPS > S-1 TOPS Deburst

4. Выделение интересующей области

Raster > Subset...

Cокращает время обработки.

5. Topo phase removal

Radar > Interferometric > Products > Topographic Phase Removal

6. Multilooking

Radar > Multilooking

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

Можно экспериментировать с величинами Number of Range Looks и Number of Azimuth Looks : 3x1, 6x2 (4x1, 8x2).

7. Goldstein filtering

Radar > Interferometric > Filtering > Goldstein Phase Filtering

8. Развертка (восстановление) фазы (Phase unwrapping)

  1. Экспорт в SNAPHU: Radar > Interferometric > Unwrapping > Snaphu Export
  2. Открыть окно терминала. Перейти в папку, куда были экспортированы файлы (cd).
  3. Открыть файл snaphu.conf
  4. Скопировать и запустить команду из snaphu.conf. Вычисления потребуют времени!
  5. Импортировать результаты в SNAP: Radar > Interferometric > Unwrapping > Snaphu import

Работа со SNAPHU

Snaphu доступен только для Linux. Устанавливается традиционно: sudo apt-get install snaphu

Для Windows понадобится использовать виртуальную машину или Cygwin.

Экспорт в SNAPHU: Target Folder — папка для экспорта данных — вводится как имя файла. Баг SNAP.

Запуск Snaphu:

  1. Открываем терминал.
  2. Переходим в папку, куда были экспортированы данные интерферограммы.
  3. Открываем в текстовом редакторе файл snaphu.conf. Копируем из него команду, указанную после строки Command to call snaphu.
  4. Вставляем команду в терминал и отправляем ее на выполнение.

Развертывание требует много памяти и времени.

В случае появления ошибки Unexpected or abnormal exit of child process следует во время экспорта в SNAPHU увеличить значения row и column с 10 (по умолчанию) до 20. Желательное решение проблемы: увеличить объем ОЗУ.

Импорт из SNAPHU:

  1. Read-Phase: загрузить снимок, полученный непосредственно перед экспортом в SNAPHU (после фильтра Голдстейна).
  2. Read-Unwrapped-Phase: выбрать hdr-файл развернутой фазы, подготовленный SNAPHU.

Ошибка импорта: couldn't read the hdr file.

Причина неизвестна. Возникает, если на шаге S1 TOPS Coregistration выбрать обе поляризации: VH и VV. Решается удалением из продукта VH-band'ов перед экспортом в SHAPHU.

Нужно проверить соответствие имени img-файла тому, что записано внутри hdr-файла и привести в соответствие. HDR и IMG должны иметь одинаковые имена.

9. Преобразование фазы в смещение

Преобразование фазы в смещение осуществляется командой SNAP: Radar > Interferometric > Products > Phase to Displacement. Результатом ее выполнения является смещение в направлении линии прямой видимости радара (Line-of-sight, LOS).

Расчет смещения в направлении линии прямой видимости радара

Расчет смещения в направлении линии прямой видимости радара. Источник

Положительные значения LOS-смещений означают движение поверхности по направлению к радару (поднятие, uplift), отрицательные — движение от радара (понижение, subsidence).

Расчет смещений командой Phase to Displacement осуществляется по формуле

$$ \Delta R = -\frac{\lambda}{4\pi}\phi_{diff} . $$

и представляет собой LOS-смещение, взятое со знаком "минус" и измеряемое в метрах.

Расчет вертикальных смещений поверхности

Расчет вертикальных смещений поверхности. Источник

Связь между смещениями имеет вид

$$ \Delta R = \Delta y \sin\theta - \Delta z \cos\theta , $$

где обозначения взяты из рисунка.

Таким образом формула для расчета вертикальных смещений имеет вид

$$ \Delta z = -\frac{\lambda}{4\pi \cos\theta}\phi_{diff} $$

Вычисление вертикальных смещений по развернутой интерферограмме включает следующие шаги:

  1. Создать новый band с формулой для вычисления смещений.
  2. Удалить из него области с низкой когерентностью.
  3. Выделить точку, которая полагается неподвижной, и узнать значение ее высоты.
  4. Вычесть значение этой высоты из всех точек изображения.

Шаги 3 и 4 выполнять не обязательно, поскольку выбор неподвижной точки представляет собой отдельную нетривиальную задачу.

1. Создание нового band'а

  1. Выбрать Band c развернутой фазой (unwrapped phase).
  2. ПКМ, Band Math.
  3. Убрать галочку Virtual (сделать слой "настоящим"). Назовем новый band: vert_disp.
  4. Нажать Edit Expression... и ввести формулу расчета вертикального смещения:
(Unw_Phase * wavelength_in_cm) / (-4 * PI * cos(rad(incident_angle)))

Unw_Phase — имя band'а с развернутой интерферограммой из текущего продукта.

incident_angle берется из Tie-Point Grids. Чтобы добавить эти значения, при вводе формулы нужно установить галочку Show Tie-Point Grid.

wavelength_in_cm — длина волны в сантиметрах. Для Sentinel-1 равна 5.55 см.

Результатом будет вертикальное смещение в см.

Введенные в band формулы можно редактировать через ПКМ, Properties > Pixel-Value Expression

Убирать галочку Virtual обязательно, т.к. виртуальные band'ы не учитываются на следующих этапах обработки (например, при геокодировании). Значок виртуального band'а помечен буквой "V" (ПКМ, Convert Band должно преобразовывать виртуальный band в реальный, но работает некорректно).

В поле Unit ввести единицы измерения, например, cm. Эти единицы будут затем отображаться на цветовой легенде в Google Earth. Добавить Unit можно и после создания band'а: ПКМ, Properties > Unit

SNAP самостоятельно выбирает какой из снимков сделать Master'ом, а какой Slave'ом на этапе S1 TOPS Coregistration. Поэтому нужно проверить временную последовательность снимков Master и Slave по метаданным. Если окажется, что Slave-снимок сделан позже Master'а, то формулу нужно откорректировать, умножив на -1.

2. Удаление областей с низкой когерентностью

В vert_disp, ПКМ, Properties, в поле valid pixel expression указываем логическое условие вида:

coh_HH_11Feb2008_04Mar2008 > 0.6

coh_HH_11Feb2008_04Mar2008 — имя band'а с когерентностью, его нужно предварительно скопировать (ПКМ, Properties, скопировать Name).

Пороговое значение когерентности варьируется и начинается примерно с 0.3.

Теперь у нас есть относительное вертикальное смещение между двумя снимками.

3. Выбор точки, смещение которой равно нулю

Предположим, что имя band'а — vert_disp, а смещение (высота) точки равно 15 (ищется по Pixel Info).

4. Вычитание значений высоты

  1. Выделить vert_disp.
  2. ПКМ, Band Math, убрать галочку Virtual. Назовем новый band: vert_disp_ref.
  3. В поле выражений записать
vert_disp - 15

Теперь у нас есть абсолютное вертикальное смещение, т.е. смещение относительно неподвижной точки.

Если интересуют только LOS-смещения, нужно попробовать удалить области с низкой когерентностью еще из развернутой фазы, т.к. сама операция Phase to Displacement не предполагает выбора параметров.

10. Геокодирование

Radar > Geometric > Terrain Correction > Range-Doppler Terrain Correction

11. Экспорт результатов

  1. В KMZ для визуализации в Google Earth: открыть для просмотра нужный band и ПКМ, Export View as Google Earth KMZ.

Перед экспортом в KMZ настроим цвет карты смещений:

  1. Во вкладке Colour Manipulation выполнить импорт палитры. Например, jet.
  2. Выполнить 95%-масштабирование (обрезает нижние 1% и верхние 4% диапазона значений отображаемой величины).

Для визуализации в QGIS удобнее использовать формат GeoTIFF: выделить конечный продукт и File > Export > GeoTIFF (или GeoTIFF/BigTIFF).

Комментарии

comments powered by Disqus