Вычислительные технологии
Том 17, № 6, 2012
Веб-сервис для построения интерферометрических
моделей ALOS L1.0
В. П. Потлпов, С.Е. Попов, А. В. Семёнов Институт вычислительных технологий СО РАН, Кемеровский филиал, Россия e-mail: [email protected], [email protected], [email protected]
Представлен веб-сервис по обработке RAW-данных радиолокационной интерферометрии уровня 1.0, построенный на базе свободно распространяемого программного обеспечения GMTSAR. Приведено описание функциональных компонентов системы, а также результатов, полученных в процессе реализации принципиальных схем двухшаговой обработки радарных снимков. Показаны примеры графического интерфейса веб-сервиса на базе технологий ASP.NET, PHP и Google API.
Ключевые слова: интерферометрия, веб-сервис, радарные снимки, фокусировка, развёртка фазы.
Введение
В последние годы для изучения природной среды всё активнее применяются радиолокационные системы с синтезируемой апертурой, главным достоинством которых является возможность получения изображений независимо от состояния атмосферы и условий естественной освещённости местности. В нашей стране для решения таких задач угледобывающей отрасли, как высокоточный учёт объёмов добычи в разрезах открытого типа, обнаружение и наблюдение просадок грунтов, методы радиолокационной интерферометрии только начинают применяться, в то время как во многих странах (США, страны Евросоюза, Австралия) мониторинг поверхности с использованием этих методов уже широко распространён. Преимущества такого вида мониторинга следующие:
— возможность построения высокоточных (в зависимости от данных до 1 см и менее) высотных моделей рельефа земной поверхности;
— всепогодность в силу высокой устойчивости радарного сигнала к атмосферным явлениям;
— регулярность вследствие размещения радарных сенсоров на космических аппаратах;
— сравнительная простота получения качественной модели рельефа и минимальные трудозатраты при обработке данных;
— относительно небольшая стоимость за 1 км2.
Большинство коммерческих продуктов, таких как ERDAS Imagime, Geomatica, ENVI, работают с фокусированными радарными данными (L 1.1). Они позволяют, манипулируя значениями растровых данных и их географической позицией, обнаруживать особенности местности, определять географические координаты объектов, извлекать линейные объекты, разрабатывать пространственные модели обработки информации (spatial modeler), переводить данные из одного формата в другой (import/export),
ортотрансформировать, составлять мозаики из изображений, получать стереоизображения и автоматически извлекать информацию по географии объектов. Несмотря на большие функциональные возможности, у коммерческих продуктов в области данных дистанционного зондирования (ДДЗ) есть один значительный недостаток — их высокая стоимость, которая в "топовых" версиях может доходить до 100 тыс. долларов. Более того, все они не позволяют работать с нефокусированными данными уровня 1.0 (L 1.0). Учитывая же стоимость одного L 1.1 снимка (нижняя граница цен начинается с 200 тыс. руб), ситуация для научно-исследовательских учреждений России становится крайне сложной. С другой стороны, необработанные ("сырые") ДДЗ уровня 1.0 (L 1.0) могут стоить на порядки дешевле. К примеру, японская служба космических данных CROSS-EX (Internet data service system dedicated to ALOS) предлагает диапазон цен за стереопару не более 40 тыс. руб. К тому же создано большое количество свободно распространяемого программного обеспечения, позволяющего вести обработку ДДЗ L 1.0 (ROI_PAC, GMTSAR, библиотеки GMT, GDAL и т.п.). Однако у таких продуктов также есть недостатки. Во-первых, их компиляции и установка в UNIX-системы за счёт подключения дополнительных библиотек и сборок являются нетривиальными. Во-вторых, ввод информации на обработку осуществляется в режиме консоли и содержит иногда огромное количество параметров, конфигураций и схем, что крайне неудобно для пользователя в плане эргономики и скорости работы. Поэтому создание на принципах WYSIWIG интуитивно понятного интерфейса для пре- и постпроцессорной обработки RAW-данных радиолокационной интерферометрии, инкапсулирующего от конечного пользователя сложное взаимодействие программных функционалов, позволяет сосредоточиться на исследуемой проблеме и является, по мнению авторов, весьма актуальной задачей.
1. Математическая модель
процесса интерферометрической обработки
Стандартный процесс интерферометрической обработки, взятый в основу и реализованный в рассмотренном веб-сервисе, состоит из нижеописанных основных этапов (рис. 1).
Фаза интерферограммы складывается из нескольких составляющих:
ф =$topo+$def+$atm+$n, (1)
где Ф^,рС) — фазовый набег за счёт обзора топографии под двумя разными углами; $def — фазовый набег за счёт смещения поверхности в период между съёмками; Ф^т — фазовый набег за счёт различия длин оптических путей в силу преломления в среде распространения сигнала; Фп — вариации фазы в результате электромагнитного шума [1, 2]. Последние две составляющие не несут информации о топографии поверхности, поэтому с помощью коррекции и фильтрации исключаются;
$topo = "Г-B—^Тz + $flat earth, (2)
Лг sin в
где Б± — перпендикулярная составляющая базовой линии, соединяющей положения космического аппарата при повторной съёмке конкретной точки земной поверхности;
Рис. 1. Этапы интерферометрической обработки
в — угол обзора поверхности при первом пролёте; Л — длина волны сканирующего излучения; r — расстояние между антенной и точкой на поверхности; z — высота поверхности над опорным эллипсоидом; Фа^ earth — предопределённая фаза, рассчитанная из модели опорного эллипсоида.
Устранение мешающей составляющей Фа^ earth необходимо для уменьшения высокочастотных переходов фазы, возникающих на интерферограмме вследствие разницы расстояния, проходимого сигналом при съёмке с разных положений, а также переходов, вносимых рельефом поверхности. Информацией для расчёта Фа^ earth может быть цифровая модель рельефа низкого разрешения, в качестве которой можно использовать генерируемые комплексом GMTSAR топографические поверхности или векторизованные топографические карты масштабов 1: 100 000, 1: 200 000. Оцифровка изолиний и отметок высот на топографических картах и дальнейшая интерполяция с целью получения однородной сетки — достаточно трудоёмкий процесс, увеличивающий время построения цифровых моделей рельефа, что также указывает в пользу создания веб-сервисов с распределением вычислительной составляющей, используя облачные вычисления. В ходе интерферометрической обработки данных разных спутников установлено, что наилучшим вариантом опорного рельефа является поверхность с постоянной высотой LDB типа ASTER1 с разрешением 30 м. В этом случае не возникает необоснованных резких всплесков фазы, а остаточные неопределённости, вызванные диапазоном регистрации фазы отражённого сигнала радаром от 0 до 2п, устраняются на этапе развёртки фа-
зы [2]. Величина смещения земной поверхности Аг, произошедшего за время между повторными съёмками, отражается в фазовой составляющей Ф^:
$def = у AR. (3)
Из уравнений (1)-(3) следует, что интерферометрическая фаза Ф содержит информацию как о рельефе, так и о смещении, при этом Фdef проявляется тем больше, чем больше значение Б±, т.е. чем дальше находится спутник при своем повторном пролёте от первого положения. Расчёт цифровых моделей рельефа на основе данных PALSAR показал, что успешное восстановление рельефа возможно при длине перпендикулярной составляющей базовой линии более 3000 м, тогда как построение карты смещений необходимо выполнять при значениях Б±, не превышающих 1000 м. В этом случае разница в углах обзора поверхности незначительна, сигнал при первом и повторном пролётах проходит по близким путям, за счёт чего дополнительно возрастает пространственная корреляция сигналов [2].
Когерентность Coh является величиной применимости пары радарных снимков для интерферометрической обработки и уменьшается в силу пространственной и временной декорреляции и шума, вносимого в сигнал при распространении и отражении:
Coh = ^4=, (4)
Vsisi ■ S2
где s1 и s2 — комплексные значения отражённого сигнала для первого (master) и второго (slave) снимков. Значение когерентности, близкое к 1, свидетельствует о стабильности интерферометрической фазы, низкая её величина -- о разрушении фазы регистрируемых сигналов [2].
2. Фокусировка
Поскольку возможность фокусировки является одним из главных преимуществ пакета СМТВАИ, по сравнению с коммерческими продуктами, отдельно стоит отметить процесс фокусировки радарных изображений, полученных сенсорами с синтезированной апертурой. Процесс фокусировки требует значительных вычислительных ресурсов, но его алгоритм довольно прост. При фокусировании БАЯ-изображения (полученного радаром с синтезированной апертурой) сначала производится последовательное суммирование упорядоченных по расстоянию сигналов вдоль следования сенсора синтетической апертуры (рис. 2). Следовательно, расстояние Я вычисляется как функция от времени Б:
Я(з) = Яо + Яо(з - во) + Я(« - зо)2 + ... + Фае! = ДЯ, (5)
2 Л
где Яо — максимальное расстояние сближения космического аппарата и точки отражателя, з0 — время максимального сближения. Параметры Я0, Я0 и Я0 необходимы для фокусировки изображения. В терминах БАЯ-обработки они называются ближними Я0. Доплеровский центроид /дс и доплеровская частота /я связаны с коэффициентами многочлена (5) как
/Ьс = —л~ /я = —Л~' (6)
Рис. 2. Геометрия БЛЯ-антенны, пролетающей над точкой-отражателем, расположенной на земной поверхности. Расстояние Я зависит от интервала времени в, измеряемого по точным координатам орбиты в координатной системе, привязанной к Земле
где Л — длина волны радара. Если ЯЛИ-изображение фокусируется с нулевым коэффициентом Доплера, то диапазон скорости в точке ¿о будет определён нулевым.
Далее производится сопоставление каждой точки на поверхности земли (широта/долгота, топография) с координатами расстояние/азимут радара. Если радиолокационное изображение фокусировалось на доплеровский центроид, отличный от нуля, то уравнение (6) может использоваться для расчёта ненулевого диапазона при этом коэффициенте Доплера:
Я = - ^. (7)
Необходимо отметить, что производная параболической аппроксимации определяет сдвиг по времени вдоль трека, который создаётся при фокусировке с ненулевым коэффициентом Доплера. Это даёт
Д* = - ^. (8)
2 Я V У
Соответствующая поправка диапазона задается в виде
ДЯ = Яд*2. (9)
Эти две поправки необходимы для перехода от географических координат к радарным в случаях, когда изображение сфокусировано с коэффициентом Доплера, отличным от нуля. В результате приведённых преобразований создаётся комплексный ЯЬС-файл, содержащий значения фазы, магнитуды и координатные привязки в системах расстояние/азимут и широта/долгота.
3. Программный функционал веб-сервиса ALOS L 1.0 processor
Веб-сервис представляет собой распределенное решение, построенное на базе технологий ASP.NET и PHP, отличительной чертой которого является распределение не только клиент-серверной части, но и сервер-уровня. Функционал сервера программно разделён на две составляющие (далее модуль ALOS Processor и модуль SAR Calc): веб-логика
и расчётная логика. Первая реализована на базе .NET Framework 4.O, развёрнута на вебсервере IIS 7.O, вторая построена на базе PHP с применением csh-интерпретатора, развёрнута на веб-сервере Apache с применением службы FTP (ОС Ubuntu Server lO.O4). В качестве расчётного ядра был выбран свободно распространяемый пакет GMTSAR, построенный на основе библиотек GMT и NetCDF, отличительной особенностью обработки которых является возможность preprocess фокусировки радарных данных ALOS уровня L l.O. GMTSAR содержит набор csh-скриптов, составляющих исполняемую среду обработки переданных на вход как пары SAR-снимков, так и пакетных картриджей из трёх и более файлов, синхронизированных по времени, в порядке возрастания временных отметок. Принципиальная схема двухшаговой обработки радарных снимков представлена на рис. 3. Описание csh-скриптов GMTSAR приведено ниже:
— align.csh — соотнесение пары SAR-снимков;
— align_batch.csh — соотнесение множества пар SAR-снимков;
Рис. 3. Принципиальная схема двухшаговой обработки радарных снимков
— baseline_table.csh — расчёт базовой линии и таблицы времени;
— cleanup.csh — очистка рабочих директорий;
— dem2topo_ra.csh — конвертация DEM в координаты расстояние/азимут;
— filter.csh — фильтрация интерферограммы, генерация изображений амплитуды, фазы и корреляции;
— fitoffset.csh — вычисление параметров аффинного преобразования;
— geocode.csh — преобразование расстояние/азимут в широта/долгота;
— grd2kml.csh — генерация kml-файла для Google Earth;
— intf.csh — генерация интерферограммы на основе одной пары SLC-изображений;
— intf_batch.csh — генерация интерферограммы на основе множества пар SLC-изоб-ражений;
— make_a_offset.csh — проведение смещения по азимуту;
— make_dem.csh — создание DEM из граней;
— pre_proc.csh — предварительная обработка (фокусировка) пары "сырых" SAR-данных;
— pre_proc_batch.csh — предварительная обработка (фокусировка) множества пар "сырых" SAR-данных;
— process2pass.csh — создание интерферограммы от начала до конца;
— process2pass_envi.csh — создание интерферограммы от начала до конца для Envisat;
— proj_ll2ra.csh — проецирование grd-файла из широта/долгота в расстояние/азимут;
— sarp.csh — фокусировка одного SAR-снимка;
— slc2amp.csh — генерация изображения амплитуды из SLC-файла;
— snaphu.csh — развёртка фазы путем использования модуля snaphu;
— update_PRM.csh — обновить значения в PRF-файле.
Каждый из перечисленных csh-скриптов принимает на входе определённое количество параметров. Например, как видно из рис. 3, скрипт pre_proc.sch оперирует с двумя файлами формата CEOS (ALOS L1.0 *.0_A), скрипт dem2topo работает с предварительно сгенерированным файлом топографии (dem.grd). Скриптам intf.csh и filter.csh необходимо указывать названия гаусс-фильтров и т. д. Логика работы всех скриптов позволяет создать единые конфигурационный файл и точку запуска, передавая только начальные параметры и сам конфигурационный файл. Для этого в пакете GMTSAR предусмотрен исполняемый скрипт process2pass.csh. Модуль SAR Calc состоит из PHP-скриптов, формирующих файл конфигурации на основе переданных клиентом параметров и запускающих на исполнение csh-скрипты GMTSAR. Для каждого пользователя веб-сервиса целесообразно создавать копию рабочей папки GMTSAR с размещением её в домашнем каталоге и производить запуск csh-скриптов от имени того же пользователя, что существенно повысит безопасность веб-сервиса и даст возможность контролировать ресурсы сервера штатными средствами UNIX (limits.conf+rlimits+groups). Более того, размещение PHP-скриптов в рабочей папке GMTSAR позволит дифференцировать виртуальные папки на веб-сервере и персонализировать настройки процесса обработки радарных данных в каждом конкретном случае. В дальнейшем авторами планируется создание отдельного веб-компонента для формирования расчётных данных. Аналогично размещению веб-ресурсов GMTSAR для каждого пользователя создаётся ftp-ресурс для выгрузки файлов L 1.0 (*.0_A) и загрузки результатов обработки. Иначе говоря, весь функционал сервер-уровня для конкретного пользователя размещён в отдельном виртуальном и физическом пространстве на сервере, что позволяет квотировать его средствами UNIX-системы (fstab+quota). Контроль и управле-
ние GMTSAR-процессами пользователей осуществляется посредством UNIX-утилиты (ps+top). Поскольку каждый пользователь запускает процесс обработки SAR-данных от своего имени, то процесс мониторинга и управления сводится к командному интерпретатору UNIX. К примеру, для принудительного завершения процесса обработки пользователя userl применяется команда killall -g < ps -f -u userl | grep process2pass.csh | cut -d ' ' -fl. Модуль SAR Calc включает следующие PHP-скрипты:
— alosl10_config.php — непосредственно создаёт файл конфигурации для csh-скрип-тов и запускает их от имени указанного в параметре пользователя; создаёт именованный UNIX-процесс;
— alosl10_rawSARdata.php — возвращает список файлов (*.0_A, RAW SAR Data), выгруженных пользователем, на FTP сервер;
— alosl10_status.php — возвращает 1 или 0, запущен ли процесс пользователя или нет соответственно;
— alosl10_log.php — возвращает log-файл процесса работы csh-скрипта process2pa-ss.csh;
— alosl10_results.php — возвращает список файлов результатов работы csh-скрипта process2pass.csh;
— alosl10_kml.php — возвращает содержимое kml-файлов, созданных csh-скриптом process2pass.csh, для отображения в Google Earth Plugin.
Рис. 4. Веб-форма для подготовки параметров обработки радарных данных ЛЬСЯ Ь 1.0 пакетом GMTSAR (фрагмент)
Модуль ALOS Processor содержит веб-форму (рис. 4) для задания параметров конфигурационного файла GMTSAR в виде интуитивно понятного GUI-интерфейса. Веб-форма разделена на секции, каждая из которых соответствует параметрам задаваемых для этапов принципиальной схемы двухшаговой обработки радарных снимков (см. рис. 3). Форма позволяет выбрать главную и зависимую радарную сцену, предварительно загруженную на FTP-сервер (файл формата RAW SAR Data, *.0_A), а также подготовить топографическую сетку согласно географическим координатам сцен. Топографическая основа генерируется при помощи online веб-сервиса GMTSAR. Функционал модуля построен на базе метода WebRequest, позволяющего получать сетевой поток ввода/вывода PHP-скриптов модуля SAR Calc.
Модуль ALOS Processor содержит веб-форму (рис. 5) для просмотра результатов работы модуля SAR Calc. Ниже представлено описание файлов результатов расчёта, созданных на базе GMTSAR с использованием csh-скрипта process2pass.csh:
— *.SLC — представляет собой SLC-изображение (Single Look Complex); четырёхбайтовые комплексные числа со знаком представлены как целые числа короткого целого без знака (двухбайтовые вещественная и мнимая части); числа масштабируются для использования полного динамического диапазона типа short integer;
— trans.dat — содержит отображение пары координат range/azimuth в координаты lon/lat топографии; топография вычислена относительно локального сферического приближения к земле, и поэтому в отображении возможны большие порядки отрицательных и положительных чисел;
Рис. 5. Веб-форма отображения результатов работы модуля ALOS Processor с возможностью представления результатов на ЭБ-картах Google Earth
— topo_ra.grd — сетка формата Netcdf со значениями расстояния, азимута и топографии;
— phase.grd — сетка формата Netcdf со значениями расстояния, азимута и модулированной интерференционный фазы (2П);
— corr.grd — сетка формата Netcdf со значениями расстояния, азимута и интерференционной когерентности (от 0 до 1);
— xphase.grd — сетка формата Netcdf со значениями расстояния, азимута и градиента фазы в направлении расстояния (радиан/пиксель);
— yphase.grd — сетка формата Netcdf со значениями расстояния, азимута и градиента фазы в направлении азимута (радиан/пиксель);
— mask.grd — сетка формата Netcdf со значениями расстояния, азимута либо NaN или 1, где NaN показывает, что данные должны быть маскированы;
— unwrap.grd — сетка формата Netcdf со значениями развернутой интерференционной фазы (радианы).
Модуль SAR Calc для каждого файла из перечня генерирует трансформированный файл в географических координатах (файлы с суффиксом "_ll", например, phase_ll.grd), файл формата postscript и PNG (расширение и .png) и файлы изображений Google Earth (расширение .kml). Возможна комбинация таких преобразований, например, файл phase_ll представляет postscript описание графических объектов фазы в географических координатах. Веб-форма модуля ALOS Processor, реализованная на базе Google API с применением JavaScript, позволяет просматривать эскизы результатов процесса интерферометрии (.png файлы) и их отображения в Google Earth (.kml файлы).
cf-xmi ver5xon="1.0u encoding="UTF-8"7>
ekml xffllns=Hhttp://earth.google.com/kml/2,1">
<oocumeitt>
<name>phase raask ll</name> Файл результата , работы SAR Cale
<CrourdOverlay>
<name>GMT Image Overlay</name> *
<Icon>
<href>http;//localhost:8082/home/user]/gintsar/csh
/intf/2999351 2919124/phase mask U .png</href=>
</Icon>
<altitudeMode>clampToGround</altitudeHode>
<LatLonBox>
<north>33.l66667</north>
<south>32.166667</south>
<east>-115.166667</east>
<west>-116,2599G9« wests-
</LatLonBox> — Область б Google Earth
<Region> . (допгатаУширота)
<LatLonAltBox>
<north>33.166667</north>
■<south>32.166667</south^
<east>-115-166667</east>
<west>-116.2599<39<:/west>
</LatLonAitBox>
</Reglon>
</GroundOverlay>
c/Document>
< t km\>_
Рис. 6. Формат генерируемого KML-файла для Google Earth
При этом модель Google Earth формируется динамически, т. е. при переходе по ссылке "View in Google Earth" полученный сетевой поток преобразуется в XML-описание (рис. 6), которое передается как строковая переменная JavaScript, использующая Google API (метод parseKml), сам же файл (.png) результата накладывается на рельеф как изображение в изометрии.
Заключение
Разработанный веб-сервис является инструментом научно-практического применения интерферометрической технологии обработки радарных изображений ALOS L 1.0, позволяющим в интерактивном online-режиме формировать расчётные процессы GMTSAR и просматривать результаты в графическом виде, что значительно снижает трудоёмкость и повышает качество анализа данных на базе данных дистанционного зондирования, сокращая тем самым время и стоимость научно-исследовательских работ.
Список литературы
[1] Филатов А.В. Обнаружение подвижек земной поверхности в зоне интенсивной нефтедобычи методами радарной интерферометрии // Вестник Югорского гос. ун-та. 2006. № 4. C. 103-109.
[2] Евтюшкин А.В., Филатов А.В. Технология построения цифровых моделей рельефа местности и оценки смещений методом радарной интерферометрии // Вестник НГУ. 2009. Т. 7, вып. 1. C. 66-72.
Поступила в 'редакцию 17 октября 2011 г., с доработки — 21 сентября 2012 г.