CC BY
127
НАУЧНОЕ ИЗДАНИЕ МГТУ ИМ. Н.Э. БАУМАНА
НАУКА и ОБРАЗОВАНИЕ
Эл № ФС 77 - 48211. Государственная регистрация №0421200025. ISSN 1994-0408
электронный научно-технический журнал
Развертка фазы радиолокационных топографических интерферограмм
# 07, июль 2012
Б01: 10.7463/0712.0423364
Шувалов Р. И.
УДК. [52-17::528.8.044.2]+519.176
Россия, МГТУ им. Н.Э. Баумана Shuvalov.R.BMSTU@mail.ru
Введение. Большое прикладное значение имеет метод построения цифровых моделей рельефа (ЦМР) поверхности Земли по данным интерферометрической съемки из космоса радиолокатором с синтезированной апертурой антенны (РСА). Это связано с возможностью получать радиолокационные изображения независимо от времени суток и погодных условий; с оперативностью получения ЦМР на интересующий район; с высокой потенциальной точностью метода. Особенно метод актуален для протяженной территории России, большую часть года скрытую облачным покровом от объективов космических оптических сенсоров.
Интерферометрический метод построения ЦМР по данным РСА заключается в проведении двух космических радиолокационных съемок интересующего участка поверхности Земли с незначительно различающимися углами наблюдения, формировании топографической интерферограммы по результатам совместной обработки полученных снимков, и извлечении из сформированной интерферограммы топографической информации [1-4].
Топографическая интерферограмма представляет собой матрицу главных (т.е. известных по модулю 2п радиан) значений разностей фаз. Для извлечения из такой интерферограммы информации о рельефе необходимо преобразовать ее в матрицу абсолютных значений фазовых разностей. Задача восстановления массива абсолютных фазовых значений по массиву главных значений фазы называется задачей развертки фазы. Развертка фазы является наиболее сложным этапом интерферометрической технологии, поскольку на получаемых интерферограммах практически всегда присутствуют разрывы, положение которых неизвестно [5-7]. Для установления положения разрывов фазы на интерферограмме необходима дополнительная информация. В космической радиолокационной топографической интерферометрии дополнительная информация включает интенсивность принятого сигнала на двух снимках, когерентность между двумя
снимками, и априорную информацию о рельефе. После включения в постановку задачи развертки фазы дополнительной информации и выбора критерия оптимальности решения задача развертки фазы становится задачей оптимизации. Полученная задача оптимизации, как правило, является нелинейной и имеет большую размерность.
Проблеме развертки фазы в космической радиолокационной топографической интерферометрии посвящено большое число зарубежных работ [4-7, 13-14], в которых описано множество методов развертки фазы, каждый из которых для выяснения положения разрывов использует лишь некоторую часть доступной дополнительной информации. Отечественных работ, рассматривающих задачу развертки фазы топографических РСА-интерферограмм, сравнительно мало, несмотря на то, что в России разрабатываются системы радиолокационного наблюдения Земли из космоса (Кондор-Э, Ресурс-Метеор 3, Аркон-2М). В работах А.И. Захарова и Л.Н. Захаровой (ИРЭ РАН) [8, 9], А.С. Леонова и Д.Д. Дарижапова (ОФП БНЦ СО РАН) [10] исследуются и сравниваются между собой различные методы развертки фазы. В работе Р.Р. Ковязина (СПбГУ ИТМО) [11] для выполнения развертки фазы интерферограммы используется метод локального интегрирования, при этом предполагается отсутствие разрывов восстанавливаемой абсолютной фазы. В работе А.В. Филатова (ЮНИИ ИТ) [12] предлагается перед разверткой фазы выполнять некогерентное накопление интерферограммы, но это снижает точность получаемого решения.
Общим подходом к снятию неоднозначности решения, связанной с неизвестным положением разрывов фазы на интерферограмме, является построение на множестве допустимых решений распределения вероятностей (степеней доверия) с учетом всей доступной информации. Вектор, компонентами которого являются две абсолютные фазовые разности, соответствующие двум взаимно перпендикулярным направлениям на плоской цифровой интерферограмме, называется градиентом абсолютной фазы. Получение распределения вероятностей на множестве положений разрывов сводится к построению распределения вероятностей градиента абсолютной фазы.
Целью настоящей работы является разработка математической модели градиента абсолютной фазы на радиолокационной топографической интерферограмме и создание метода развертки фазы, опирающегося на эту модель, и учитывающего всю доступную дополнительную информацию.
Постановка задачи. Обобщенная постановка задачи развертки фазы на плоскости имеет следующий вид:
ф[ ^ ( х, У ), а ( х, у)]^ шт, (1)
^х[ ? ( х, у ) + а ( х, у )] = 0,
?(X,у) = Ж[Ур(х,у)] ,
V/ х, у ) = 5 (х, у) + а (х, у), (X, у )єОс Я2, где Ф[] - подходящий регуляризирующий функционал; V - набла-оператор; Ж [•] - оператор свертки по модулю 2п радиан; (р(х,у) - заданное скалярное поле главного значения фазы; /( х, у) - искомое скалярное поле абсолютной фазы; а (х, у) - неизвестное добавочное
векторное поле; О - замкнутое связное ограниченное множество на декартовой плоскости.
Если интерферограмма представлена в цифровом виде, то постановка (1) может быть сформулирована в терминах теории транспортных сетей. Интерферограмме ставится в соответствие связный ориентированный граф О = (V, Е) (здесь V - множество вершин, Е -
множество дуг), являющийся конечной целочисленной решеткой относительно декартовой системы координат на плоскости (рис. 1).
а б
Рис. 1. Матрица главного значения фазы (интерферограмма) (а) и ассоциированный с ней
ориентированный граф (б).
При этом каждые четыре попарно-смежных пикселя интерферограммы соответствуют некоторой вершине графа, а каждая пара смежных пикселей соответствует паре противоположно ориентированных дуг. Каждой вершине I е V приписывается интенсивность е К, равная величине фазового остатка, вычисленного для
соответствующей четверки попарно-смежных пикселей интерферограммы Ф = {<рт п} :
ц (т, п) = 8Х (т, п) + 8Г (т, п +1) - 8Х (т +1, п) - 8Г (т, п),
5Т (т „) = Ж [рт+1,„ - рт „ ] , 5Х (т „) = Ж [^,„+1 - Рт„]
Если д > 0, вершина называется источником; если д < 0, вершина называется стоком; если д = 0, вершина называется нейтральной. Пропускная способность каждой дуги
пропускными способностями дуг и с определенными на дугах функциями стоимости называется транспортной сетью. Сетевая модель рассматривает распределение в сети потока некоторой субстанции, перетекающей по дугам из вершин-источников через нейтральные вершины в вершины-стоки.
Задача развертки фазы в сетевой постановке есть задача поиска потока минимальной стоимости
где д - величина потока по дуге (/', ])е Е. Функции стоимости с^ (•) нуждаются в
предварительном построении по данным измерений и априорной информации о решении. В различных приложениях развертки фазы функции стоимости могут определяться по-разному. В общем случае функции стоимости для задачи развертки фазы в сетевой постановке могут быть получены на основе распределения вероятностей абсолютной фазовой разности:
параметров; А - абсолютная фазовая разность; Р (-п < А < п) - вероятность непрерывности
Ниже установлено параметрическое распределение вероятностей абсолютной фазовой разности, позволяющее вычислять по формуле (3) функции стоимости для задачи развертки фазы радиолокационных топографических интерферограмм. Поток минимальной стоимости, являющийся решением задачи развертки фазы в сетевой постановке (2)-(3), соответствует положению разрывов фазы на интерферограмме, наиболее согласующемуся с имеющимися данными.
предполагается неограниченной. Пусть каждой дуге (і, у)є Е поставлена в соответствие скалярная функция с у (д), определяющая стоимость протекания потока величины д є Я, д > 0 по этой дуге. Граф О = (V, Е) с заданными интенсивностями вершин, с заданными
(2)
(і, У )єЕ
Е %- Е ь =^ viєV,
{у. (і,У)єЕ} {у. (у,і)єЕ}
, к є 7 , (і, у)є Е, (3)
где к - величина потока по дуге (кратность разрыва фазы); 6у
вектор значений
фазы; Р (—п + 2пк < А < п + 2пк) - вероятность наличия разрыва кратности к.
Математическая модель. Наблюдаемые значения интерферометрической фазы, когерентности и интенсивности зависят от большого числа факторов, в том числе от случайных факторов. Поэтому в настоящей работе при математическом моделировании используется аппарат теории вероятностей и математической статистики. Наблюдаемые значения трактуются как реализации случайных величин с известными законами распределения. На основе байесовского подхода в настоящей работе разработана математическая модель градиента абсолютной фазы на радиолокационной топографической интерферограмме. Байесовский подход позволяет использовать данные измерений (главное значение фазы, когерентность, интенсивность) совместно с априорной информацией (статистическими характеристиками рельефа покрытой съемкой местности) (рис. 2, рис. 3). Разработанная модель представляет собой пару параметрических распределений вероятностей абсолютных фазовых разностей:
ад
^ PN ( Ах — АТх , р) Р (8х I АТх , р) ? (I I АТх ) ? (АТх ) dАTX
p (Ах 18Х, р, I) = ^-----ад------------------------------------------------, (4)
^ Р (8х 1 Атх,р) р (1 1 Атх, Ату ) р (АТХ ) dАтх
^ рм ( Ау Ату ,р) р (8у 1 Ату,р) р ( Ату ) dАту
Р (М5Г, р)=-ад-------ад----------------------------------,
^ р (8У 1 АТУ , р) ? (АТУ ) dАTY
-ад
где р (•) - плотность распределения вероятностей абсолютной фазовой разности; Ах, АY -абсолютные фазовые разности по направлениям наклонной дальности (X) и азимута (У), характеризующие локальный наклон фазового рельефа; 8х ,8У - относительные фазовые разности; I - интенсивность принятого радиолокационного сигнала; р - когерентность; Атх, Ату - физические (т.е. полезные, не искаженные шумом) фазовые разности. Предложенная модель (4) состоит из нескольких компонентов:
1) р¥! (•) - плотность распределения вероятностей фазового шума;
2) р (81 Ат , р) - функция правдоподобия физической фазовой разности по наблюдаемой относительной фазовой разности 8;
3) р (11 Атх) - функция правдоподобия физической фазовой разности Атх по наблюдаемой интенсивности I радиолокационного сигнала;
4) р (Атх, Ату ) - плотность априорного совместного распределения вероятностей физических фазовых разностей Атх, Ату .
Предложенное априорное совместное распределение вероятностей физических фазовых разностей (рис. 2 справа, рис. 3, табл. 1)
[Рт (£х (Атх ) (А тх, Ату )) ^ (Атх, Ату ), Атх < Атх
PA (ATX , ATY )
0, Atx — Atx
включает в себя преобразование, связывающее локальные наклоны (производные по пространственным координатам) рельефа подстилающей поверхности с локальными наклонами фазового рельефа на интерферограмме, полученное из уравнений геометрии съемки:
(A )= Яг0sin2 (Yo )ATX__________ (5)
x TX 4^B11 Ars + Яг0 sin ( Yo ) cos ( Yo ) ATX ’
g (д A )= Яг0 Ars sin (Yo )Aty
gY \ ATX , ATY /
4^B11 ArSAa + Яг0Aa sin ( y0 ) cos ( Yo) ATX ’
А /а * \ A I \ A* 4n| B1 Ars
atx g(atx; +ад), aty є(-ад; +ад), atx = --
(го)со§ (/о)’
и априорное совместное распределение вероятностей локальных наклонов рельефа подстилающей поверхности, полученное путем функциональной аппроксимации экспериментальных гистограмм (рис. 2 слева):
Рт (£х,) = 3395ехр(-4^х + g2r )14) , gх = tan(ах), = ^п(ау ) .
Таблица 1.
Значения параметров, принятые при построении графиков на рис. 2-3.
Параметр Значение
Длина волны РСА, Я 0,000057 км
Наклонная дальность до центра кадра, г0 1027 км
Угол наблюдения для центра кадра, у0 40°
Перпендикулярная компонента базовой линии, Б± 0,109 км
Размер пикселя по направлению наклонной дальности, Аг3 0,023 км
Размер пикселя по направлению азимута, Аа 0,021 км
Коэффициент некогерентного накопления В
Рис. 2. Априорное распределение топографического градиента (слева) и совместное априорное распределение физических фазовых разностей (справа).
Р 0,045 0,04 0,035 0,03 0,025 0,02 0,015 0,01 0,005 0
V
1
!
X И и
\1
7
-Зя
-2л
2л
Ъл
Аг, рад
Рис. 3. Априорное распределение вероятностей физической фазовой разности Ат по направлению наклонной дальности (X) и по направлению азимута (У).
Разработанная функция правдоподобия физической фазовой разности Атх по интенсивности радиолокационного сигнала (рис. 4-5):
р1 (11 Атх )
1
Ь
Г( Ь )1 М (Атх )
1Ь 1 exp
- Ь
I
м (Атх)
I > 0, Ь > 1
включает в себя математическое ожидание наблюдаемой интенсивности:
M (А^ )= M (Атх (а ))= M1 (aх\ат = 0), MI (aх ,ат ) = CІG (^ а) + D, учет эффектов переналожения и радиолокационной тени:
I (у-П 2,а), «х ^у-П 2 4(ах,ат) = ^(ах,ат), у-П2<ах <У,
I (у,ат ), У^ах
и выражение для наблюдаемой средней интенсивности сигнала:
I (ах, ат ) = (ах, ат ) а (ах, ат ) ,
SF (ах,ат) = АaАr^ 1 (ах -у + ж/2)2 + “^п2 (у)а^ +1^,
(в) = У (в)(wc ехр (-уи2в2 ) + wi (1 + в2 ) Р + wd ехр (-в') СОБ0Л (в)),
0.2w (1 - w) (1 - w)2 2 , ч , ч2 г ,
----- -----, wd = --—, W = w + 0.2w (1 - w ) + (1 - w ), w е[0;1].
а0
wc =----------, wi =---------------------------, wd =
c ^ i W W
в (ах ,ат ) =
агссоБ
^ tan(ах)б1п(у) + соб(у) ^
лУ1ап2(а7)+"ьап^(аа)+Т
где ах ,аг - углы наклона рельефа по направлениям дальности и азимута; J(Атх, Атг ) -якобиан преобразования (5); B1 - перпендикулярная составляющая базовой линии; I -наблюдаемое значение интенсивности радиолокационного сигнала; Аг5 , Аa - размеры пикселя радиолокационного снимка по направлениям наклонной дальности и азимута; Л -рабочая длина волны радиолокатора; у0, г0 - угол наблюдения и наклонная дальность,
соответствующие центру кадра; Ь - количество независимых наблюдений; Г() - гамма-функция Эйлера; MI (ах ,аг) - предсказываемое радиометрической моделью значение средней интенсивности; SF (ах ,аг) - площадь ячейки на подстилающей поверхности, соответствующей данному пикселю на снимке; а0 (в) - удельная эффективная площадь рассеяния; в - угол падения радиолокационного сигнала; C, D, и, Р, w - параметры разработанной радиометрической модели; у - угол наблюдения при съемке.
Еще одним компонентом предложенной функции правдоподобия физической фазовой разности является функция правдоподобия p (5 \ Ат, р) физической фазовой разности Ат по относительной фазовой разности 5. Плотность распределения вероятностей относительной
фазовой разности р (8 | Лт, р) зависит лишь от главного значения разности 8 -Лт когерентности р . Эта зависимость представлена графически на рис. 6.
I 0,01 0,009 0,008 0,007 0,006 0,005 0,004 0,003 0,002 0,001 0
у=о
7=0.024
-ч II сэ Л 10 \ ^ II '--ч
0,98 3,08 5,18 7,28 9,38 11,48 13,58 Д^, рад
Рис. 4. Семейство функций правдоподобия Ь (Лтх ) = р} (11 Лтх ), полученное варьированием
наблюдаемой интенсивности I.
Г 15 "
Рис. 5. Зависимость функции правдоподобия физической фазовой разности от наблюдаемой
интенсивности радиолокационного сигнала.
л
с
го
тз
ш
а*
о
с
2?
«Ё
Ь
Рис. 6. Плотность распределения вероятностей главного значения разности 8 -Лт как функция когерентности. Количество независимых наблюдений Ь = 2 .
0.95 0.9 0.Й5 0.8 0.75 0.7 0.65 0.6 0.55 0.5
О 0.2 0.4 0.6 О.в 1
Р
Рис. 7. Вероятность непрерывности фазы по направлению азимута Р (-ж < Лу кж^у, р) как функция относительной фазовой разности 8у и когерентности р .
1
Рис. 8. Вероятность непрерывности фазы по направлению наклонной дальности Р (-п < Ах <п\8х, р, I) как функция относительной фазовой разности 5Х, когерентности р
и интенсивности I.
Разработанная модель (4) позволяет включить доступную дополнительную информацию в постановку задачи (2) через функции стоимости (3). Модель не привязана к какому-либо методу развертки фазы и имеет самостоятельную ценность, т.к. позволяет оценивать вероятности разрывов фазы различной кратности на топографической РСА-интерферограмме по имеющимся данным измерений и априорной информации о рельефе местности (рис. 7-8).
Разработанное распределение вероятностей локального наклона фазового рельефа на радиолокационной топографической интерферограмме (4), по сравнению с известными ранее результатами [6, 7, 13, 14], наиболее полно учитывает имеющуюся информацию. Предложенная модель обобщает модель работы [6] по трем направлениям: учет физических разрывов фазы (т.е. разрывов, обусловленных топографией и геометрией съемки), учет интенсивности принятого радиолокационного сигнала, учет априорного распределения вероятностей топографического градиента. Разработанная модель (4) включает ряд разработанных «узкоспециализированных» моделей: модель априорной информации,
радиометрическую модель, модель формирования интерферограммы и модель фазового
шума. Каждая из этих моделей-компонентов может в дальнейшем дорабатываться независимо от других. Предложенная модель (4) представляет собой удобный инструмент интерпретации интерферограммы по доступной дополнительной информации и может использоваться совместно с различными алгоритмами развертки фазы в космической радиолокационной топографической интерферометрии.
Метод развертки фазы. В настоящей работе разработан метод решения задачи развертки фазы в постановке (2) с выпуклыми неотрицательными функциями стоимости, названный методом независимых диполей. Метод независимых диполей представляет собой модификацию известного в теории транспортных сетей алгоритма последовательного поиска кратчайших путей (англ. Successive Shortest Path Algorithm) [15] и заключается в последовательном выделении диполей (т.е. пар «источник-сток») при помощи поиска путей минимальной стоимости и пропускании вдоль найденных путей потоков единичной величины. Решение задачи развертки фазы в постановке (2) с функциями стоимости (3) эквивалентно реконструкции наиболее вероятной в смысле распределения (4) системы разрывов фазы, имеющихся на интерферограмме. Разработанный алгоритм отличается от алгоритма последовательного поиска кратчайших путей следующими особенностями: 1) пропускание отрицательных потоков наряду с положительными; 2) построение искомого потока в два этапа: на первом этапе пропускаются потоки длиной менее заданной величины р, а на втором - все оставшиеся; 3) использование предположения, согласно которому искомый поток минимальной стоимости представляет собой совокупность потоков единичной величины, длина каждого из которых существенно меньше линейного размера сети; 4) использование решения вспомогательной задачи небольшой размерности, полученной пространственным сжатием данных исходной задачи, для приближенного поиска пространственно протяженных потоков.
Вычислительный эксперимент. Точность метода развертки фазы оценивалась экспериментально путем сравнения результата работы алгоритма этого метода с эталонным результатом. Результат работы алгоритма в нашем случае представляет собой матрицу. Использовались следующие характеристики уклонения матрицы данных от эталонной матрицы: 1) среднее (по множеству элементов матрицы) уклонение; 2) средний модуль уклонения; 3) максимальное уклонение. Для более точного описания характера различия между матрицами использовались: 1) матрица пространственного распределения уклонения;
2) гистограмма уклонения. Поскольку развертка фазы является одним из этапов интерферометрической технологии построения ЦМР, ошибка развертки фазы входит составной частью в ошибку результирующей ЦМР. Поэтому, различные методы развертки фазы можно сравнивать по точности результирующей ЦМР. Этот подход к сравнительному
анализу точности методов развертки фазы с точки зрения приложений является наиболее естественным. Анализ точности проводился на двух сюжетах.
Интерферометрическая пара первого сюжета (рис. 9) получена по результатам съемки РСА БЯ8-1/2 в начале 1993 года и покрывает небольшой участок национального парка «Долина смерти» на западе США. Участок представляет собой засушливый район с редкой растительностью. Размеры участка составляют приблизительно 9х9 км. Размеры каждого из снимков пары составляют 1138х2326 пикселей. Значения основных параметров представлены в таблице ниже (табл. 3).
Таблица 3.
Значения параметров формирования интерферограммы.
Параметр Аг, м Аа, м Я, м Б±, м г0, км
Значение 7,904 3,981 0,057 133,392 857,680
В таблице 3 приняты обозначения: Аг, Аа - размер пикселя по направлениям дальности и азимута (по горизонтали и вертикали на рис. 9) соответственно; Я - рабочая длина волны РСА; Б1 - абсолютная величина перпендикулярной компоненты базовой линии; г0 -
наклонная дальность до центра кадра.
Снимки интерферометрической пары были пространственно совмещены, после чего была выделена прямоугольная область перекрытия. По области перекрытия были вычислены топографическая интерферограмма (рис. 10), матрица когерентности (рис. 11) и матрица интенсивности (рис. 12). Для снижения интенсивности шума были выполнены некогерентное накопление с коэффициентом 5 по направлению азимута и пространственная фильтрация. Затем была вычислена матрица пространственного распределения сингулярных точек (рис. 13). На рис. 14 приведено псевдоцветное изображение, полученное в модели «ЯСВ», объединяющее в себе информацию, содержащуюся в интерферограмме, в матрице когерентности и в матрице интенсивности. Изображение (рис. 14) показывает, что матрицы когерентности и интенсивности (рис. 11-12) содержат информацию о пространственном положении разрывов фазы на интерферограмме (рис. 10). На псевдоцветном изображении (рис. 14) разрывы фазы имеют вид «языков пламени» (оранжевый цвет). Разработанная в настоящей работе математическая модель (4) позволяет извлекать эту информацию и переводить ее в удобную для практического использования форму. Исходными данными для вычисления вероятности разрыва фазы заданной кратности, согласно разработанной
модели(4), являются наблюдаемая фазовая разность 8, когерентность р и интенсивность I (рис. 10-12).
Рис. 9. Интерферометрическая пара снимков (амплитудная составляющая).
Поскольку вероятности требуется вычислять для каждой пары смежных пикселей интерферограммы, а интерферограмма может иметь очень большие размеры, целесообразно предварительно построить вспомогательные таблицы: вычислить искомые вероятности в узлах сетки пространства параметров (рис. 7-8). Путем интерполяции по полученным таблицам (рис. 7-8), строятся матрицы пространственного распределения вероятностей разрывов и вероятности непрерывности фазы (рис. 15-20). Далее по полученным распределениям вероятностей в соответствии с формулой (3) были построены функции
стоимости и осуществлен переход к задаче поиска потока минимальной стоимости (2). Найденный предложенным методом независимых диполей, поток минимальной стоимости, позволил восстановить неизвестные матрицы абсолютных фазовых разностей, по известным матрицам относительных фазовых разностей (рис. 21-22). Матрицы абсолютных фазовых разностей дали матрицу абсолютной (развернутой) фазы (рис. 23).
я
.0
-0.5
1оо
Рис. 11 . Матрица когерентности.
Рис. 10. Топографическая интерферограмма.
№..
У*
Рис. 12. Матрица интенсивности радиолокационного сигнала.
Рис. 13. Пространственное распределение сингулярных точек. Изображение имеет три градации яркости: черный цвет - сингулярная точка интенсивности -1; серый цвет -отсутствие сингулярной точки; белый цвет - сингулярная точка интенсивности +1.
Рис. 14. Псевдоцветное изображение, полученное в модели «RGB»: интенсивность радиолокационного сигнала ^ красный канал; главное значение интерферометрической фазы ^ зеленый канал; когерентность ^ синий канал.
Далее на основе полученной матрицы абсолютной фазы (рис. 23) с учетом известных значений параметров съемки была построена цифровая модель рельефа (рис. 24-25). Построенная ЦМР (рис. 24-25) сравнивалась с эталонной ЦМР («USGS NED 30 meter DEM»). По результатам сравнения была построена матрица пространственного распределения погрешности полученной ЦМР (рис. 26).
Рис. 15. Пространственное распределение вероятности отрицательного разрыва фазы
(кратности -1) по направлению азимута.
0.98
0.77
0.56
Рис. 16. Пространственное распределение вероятности непрерывности фазы по направлению
азимута.
Рис. 17. Пространственное распределение вероятности положительного разрыва фазы
(кратности +1) по направлению азимута.
Рис. 18. Пространственное распределение вероятности отрицательного разрыва фазы (кратности -1) по направлению наклонной дальности.
Рис. 19. Пространственное распределение вероятности непрерывности фазы по направлению
наклонной дальности.
Рис. 20. Пространственное распределение вероятности положительного разрыва фазы (кратности +1) по направлению наклонной дальности.
Рис. 21 . Матрица относительных фазовых разностей по направлению наклонной дальности (наверху), матрица абсолютных фазовых разностей (в центре), и их небольшие фрагменты
(внизу).
Рис. 22. Матрица относительных фазовых разностей по направлению азимута (наверху), матрица абсолютных фазовых разностей (в центре), и их небольшие фрагменты (внизу).
Рис. 23. Матрица абсолютной фазы (значения выражены в радианах), полученная путем развертки фазы методом независимых диполей.
Рис. 24. Трехмерное представление полученной цифровой модели рельефа.
167 м 1570 м
Рис. 25. Цифровая модель рельефа, полученная с участием метода независимых диполей.
Рис. 26. Матрица пространственного распределения погрешностей ЦМР, полученной с
участием метода независимых диполей.
Был проведен сравнительный анализ точности различных методов развертки фазы (табл. 4, рис. 27).
Таблица 4.
Погрешности ЦМР по первому сюжету.
Алгоритм бты, м бМАХ , м М [б], м ,м а, м
МНКБПФ -426 419 -39,533 90,856 108,740
МНК СГ -423 418 -39,616 91,117 108,940
ВМНК ИП -328 431 -31,085 81,270 97,145
МФГ -332 443 -40,161 77,875 90,707
МРП -213 195 14,914 32,325 38,787
БКАРНи -148 164 -5,599 29,263 34,618
МНД -127 168 -5,135 29,531 35,095
В табл. 4 приняты следующие обозначения: МНК БПФ - метод наименьших квадратов без взвешивания, реализованный на основе быстрого преобразования Фурье; МНК СГ - метод наименьших квадратов без взвешивания, реализованный на основе итерационного метода сопряженных градиентов; ВМНК ИП - метод наименьших квадратов, использующий в качестве весовых коэффициентов значения когерентности, реализованный на основе итерационного метода Пикарда; МФГ - метод функций Грина; МРП - метод растущих пикселей; БКАРНи - итерационный метод поиска потока минимальной стоимости, использующий матрицу интенсивности и матрицу когерентности [7]; МНД - метод независимых диполей, разработанный в настоящей работе; бшм - минимальное уклонение
экспериментальной матрицы от эталонной матрицы; бМ^ - максимальное уклонение
экспериментальной матрицы от эталонной матрицы; М [б] - среднее уклонение
экспериментальной матрицы от эталонной матрицы; М [|б|] - средний модуль уклонения
экспериментальной матрицы от эталонной матрицы; а = среднеквадратическое отклонение уклонения от среднего значения.
М
с
н
о
н
я
а-
Ё
и
Н
О
0,11
0,1
0,09
0,08
0,07
0,06
0,05
0,04
0,03
0,02
0,01
0
1
1 \
/? N й \
К)
/ '■^1 ““ \
/ |/ V
т !
\ \
V
-400 -350 -300 -250 -200 -150 -100 -50 0 50 100 150 200 250 300 350 400
погрешность, м.
■МНКБПФ----ВМНКИП ----МНКСГ МРП--------ЗЫАРНи МФГ-------МИД
Рис. 27. Относительные частоты погрешностей ЦМР, полученных различными методами.
Интерферометрическая пара второго сюжета получена по результатам съемки РСА КаёагеаМ в начале 2000 года и покрывает участок национального парка «Долина смерти» на западе США. Участок представляет собой засушливый район с редкой растительностью. Размеры участка составляют приблизительно 71х107 км. Размеры основного снимка пары составляют 6129х20608 пикселей. Значения основных параметров съемки представлены в таблице ниже (табл. 5).
Значения параметров формирования интерферограммы.
Таблица 5.
Параметр Ат, м Аа, м X, м Б1, м г0, км
Значение 11,596 5,190 0,057 108,960 992,480
Снимки интерферометрической пары были пространственно совмещены, после чего была выделена прямоугольная область перекрытия. По области перекрытия были вычислены топографическая интерферограмма (рис. 28), матрица когерентности (рис. 29) и матрица интенсивности. Для снижения интенсивности шума были выполнены некогерентное
накопление с коэффициентом 2 по направлению азимута и пространственная фильтрация. Полученная интерферограмма (рис. 2B) имеет области полной декорреляции (области низкой когерентности на рис. 29). Число сингулярных точек на полученной интерферограмме равно 188276 (рис. 30). Mатрица абсолютной фазы, построенная методом независимых диполей, представлена на рис. 31, а полученная на ее основе ЦMP, показана на рис. 32-34. При оценивании погрешностей полученная ЦMP сравнивалась с эталонной ЦMP («USGS NED 30 meter DEM»). На рис. 35 показано пространственное распределение погрешности полученной ЦMP. Основные числовые характеристики точности полученной ЦMP приведены в табл. 6. На рис. 36 построен график относительной частоты погрешности полученной ЦMP. Сравнение проводилось лишь с методом наименьших квадратов (табл. 6, рис. 36). Mетод растущих пикселей в силу наличия областей полной декорреляции не дает адекватного результата вовсе. Алгоритм SNAPHU не отрабатывает из-за нехватки вычислительных ресурсов. Mетод независимых диполей отработал примерно за 50 минут на персональной ЭВM: Intel(R) Core(TM)i7 CPU 3.33 ГГц, RAM 8Гб.
Таблица 6.
Погрешности ЦЫР по второму сюжету.
Алгоритм 5MIN , м 5MAX , м M [5], м M 5 ,м а, м
MHK БПФ -15B1 BB3 -96 231 2B0
MHA -919 1160 9B 1BB 21B
В ходе вычислительного эксперимента на интерферограммах небольшого размера (менее 2000х2000 пикселей) предложенный метод независимых диполей не уступил по точности лучшему из известных нам алгоритмов «БКАРНи». Алгоритм «БКАРНи» практически не применим при обработке больших интерферограмм в силу высоких требований к вычислительным ресурсам. По сравнению же с методом наименьших квадратов, также работающим на больших интерферограммах, метод независимых диполей оказался более точным.
Разработанный метод позволяет сокращать время обработки за счёт возможного снижения точности получаемого решения. Это важно при обработке больших интерферограмм. При этом в получаемом решении могут возникать погрешности в виде линейных разрывов.
Рис. 28. Интерферограмма размером 6116х10154 пикселей.
Рис. 30. Пространственное распределение сингулярных точек. Общее число сингулярных
точек 188276.
Рис. 32. ЦМР, полученная с участием метода независимых диполей.
Рис. 33. ЦМР, полученная с участием метода независимых диполей, в трехмерном
представлении.
Рис. 34. Участок, полученной методом независимых диполей ЦМР, в трехмерном
представлении.
Рис. 35. Матрица пространственного распределения погрешностей полученной ЦМР.
Г олубые участки не связаны с методом развертки фазы, и скорее всего обусловлены либо атмосферными эффектами, либо другими источниками погрешности интерферометрической ЦМР. Локальные синие пятна связаны либо с погрешностями метода развертки фазы, либо с декорреляцией снимков интерферометрической пары.
0,14 0,12
н 0,1
О
£ о.оз
5 о,об
* 0,04
0,02
0 -1
Рис. 36. Относительные частоты погрешностей полученной ЦМР.
Заключение. Задача развертки фазы на плоскости является сложной математической проблемой, имеет большое прикладное значение и до настоящего времени активно исследуется. Сложность заключается в неоднозначности решения, обусловленной неизвестным положением на интерферограмме разрывов фазы, и большой размерности задачи. В настоящей работе на основе байесовского подхода разработана математическая модель градиента абсолютной фазы на радиолокационной топографической интерферограмме. Разработанная модель представляет собой параметрическое распределение вероятностей локального наклона фазового рельефа и позволяет включить доступную дополнительную информацию в постановку задачи развертки фазы для снятия неоднозначности, связанной с неизвестным положением разрывов. Модель обобщает известные ранее результаты по трем направлениям: учет разрывов фазы, вызванных рельефом местности и геометрией съемки; учет интенсивности принятого радиолокационного сигнала; учет статистических характеристик рельефа. Предложен метод развертки фазы интерферограмм. Метод является модификацией известного в теории транспортных сетей алгоритма последовательного поиска кратчайших путей и решает задачу развертки фазы в сетевой постановке, с использованием функций стоимости дуг, построенных на основе разработанной модели градиента абсолютной фазы. Проведенный вычислительный эксперимент показал высокую точность предложенного метода и его способность решать задачи большой размерности. Полученные результаты имеют большое
погрешность, м.
значение для построения цифровых моделей рельефа Земли по данным интерферометрических радиолокационных измерений из космоса.
СПИСОК ЛИТЕРАТУРЫ
1. Graham L.C. Synthetic interferometric radar for topographic mapping // Proc. IEEE, vol. 62, pp. 763-768, 1974.
2. Zebker H.A., Goldstein R.M. Topographic mapping from interferometric SAR observations // J. Geophys. Res., vol. 91, pp. 4993-4999, 1986.
3. Bamler R. Digital terrain models from radar interferometry // Photogrammetric week ‘97, pp. 93105, Wichmann Verlag, Heidelberg, 1997.
4. Rosen P. et al. Synthetic aperture radar interferometry // Proceedings of the IEEE, vol. 88, no. 3, pp. 333-382, March 2000.
5. Goldstein R.M., Zebker H.A., Werner C.L. Satellite radar interferometry: Two-dimensional phase unwrapping // Radio Sci., vol. 23, no. 4, pp. 713-720, 1988.
6. Carballo G.F., Fieguth P.W. Probabilistic cost functions for network flow phase unwrapping // IEEE Transactions on Geoscience and Remote Sensing, vol. 38, no. 5, pp. 2192-2201, 2000.
7. Chen C.W. Statistical-cost network-flow approaches to two-dimensional phase unwrapping for radar interferometry: PhD thesis, Stanford University, 2001.
8. Zakharov A.I. On the way of estimation of the reliability of the interferometric pixels for correct phase unwrapping in the DEM generation // Proceedings of the FRINGE‘99 Conference, pp. 1-8, Liege, Belgium, 1999.
9. Zakharova L.N. Comparison of global and local approach to phase unwrapping for rugged terrain // Proceedings of the FRINGE 2003 Workshop, pp. 1-5, Frascati, Italy, 2003.
10. Леонов А.С., Дарижапов Д.Д. Исследование методов развертки фазы для интерферометрической обработки радиолокационных данных // Современные проблемы дистанционного зондирования Земли из космоса.: Тез. докл. Всерос. конф. - Москва, 2005. -С. 31.
11. Ковязин Р.Р. Двумерное восстановление фазы интерферограмм // Проблемы когерентной и нелинейной оптики, СПб, 2000. - С. 267-275.
12. Филатов А.В. Метод обработки комплексных радиолокационных интерферограмм в условиях высокой временной декорреляции: Дис. канд. физ.-мат. наук. - Барнаул, 2009.
13. El-taweel G.S. Enhanced model for generating 3-D images from RADAR interferometric satellite images // WSEAS Transactions on Environment and Development, vol. 3, issue 3, pp. 5964, 2007.
14. Refice A. et al. Weights determination for minimum cost flow InSAR phase unwrapping // Proceedings of the IGARSS‘99 Conference, vol. 2, pp.1342-1344, Hamburg, Germany, 1999.
15. Ahuja R.K, Magnanti T.L., Orlin J.B. Network flows: theory, algorithms and applications. Prentice Hall, 1993.
SCIENTIFIC PERIODICAL OF THE BAUMAN MSTU
SCIENCE and EDUCATION
EL № FS 77 - 4821 1. №0421200025. ISSN 1994-0408
electronic scientific and technical journal
Phase unwrapping of radar topographic interferograms
# 07, July 2012
DOI: 10.7463/0712.0423364
Shuvalov R.I.
Russia, Bauman Moscow State Technical University
Shuvalov.R.BMSTU@mail.ru
Two-dimensional phase unwrapping problem was considered. Phase unwrapping is a key step of generating digital elevation model with the use of topographic radar interferometry. Difficulty of phase unwrapping problem is caused by unknown location of discontinuities on the interferogram and large interferogram sizes. On the base of Bayesian approach the mathematical model of the phase slopes on radar topographic interferograms was developed. Phase unwrapping may be considered as an optimization problem. For this statement of the problem a mathematical method based on the developed mathematical model was proposed. The computational experiment confirmed practical value of the obtained results.
Publications with keywords: two-dimensional phase unwrapping, topographic interferogram, topographic SAR interferometry, digital elevation model (DEM), synthetic aperture radar (SAR)
Publications with words: two-dimensional phase unwrapping, topographic
interferogram, topographic SAR interferometry, digital elevation model (DEM), synthetic aperture
radar (SAR)
References
1. Graham L.C. Synthetic interferometer radar for topographic mapping. Proceedings of the IEEE, 1974, vol. 62, no. 6, pp. 763-768. DOI: 10.1109/PROC.1974.9516
2. Zebker H.A., Goldstein R.M. Topographic mapping from interferometric synthetic aperture radar observations. Int. Geoscience and Remote Sensing Symp. (IGARSSS5) Digest. 1985, pp. 113-117.
3. Bamler R. Digital terrain models from radar interferometry. Photogrammetric Week ’9?. Heidelberg, Wichmann Verlag, 1997, pp. 93-105.
4. Rosen P.A. Synthetic aperture radar interferometry. Proceedings of the IEEE, 2000, vol. 88, no.
3, pp. 333-380.
5. Goldstein R.M., Zebker H.A., Werner C.L. Satellite radar interferometry: Two-dimensional phase unwrapping. Radio Science, 1988, vol. 23, no. 4, pp. 713-720.
6. Carballo G.F., Fieguth P.W. Probabilistic cost functions for network flow phase unwrapping. IEEE Transactions on Geoscience and Remote Sensing, 2000, vol. 38, no. 5, pp. 21922201. DOI: 10.1109/36.868877.
7. Chen C.W. Statistical-cost network-flow approaches to two-dimensional phase unwrapping for radar interferometry. PhD Thesis. Stanford Univ., 2001.
8. Zakharov A.I. On the way of estimation of the reliability of the interferometric pixels for correct phase unwrapping in the DEM generation. Proc. FRINGE‘99 Conf. Liege, Belgium, 1999, pp. 1-8.
9. Zakharova L.N. Comparison of global and local approach to phase unwrapping for rugged terrain. Proc. FRINGE 2003 Workshop. Frascati, Italy, 2003. pp. 1-5.
10. Leonov A.S., Darizhapov D.D. Issledovanie metodov razvertki fazy dlia interferometricheskoi obrabotki radiolokatsionnykh dannykh [Research of the methods of phase scanning for interferometric processing of radar data]. Sovremennyeproblemy distantsionnogo zondirovaniia Zemli iz kosmosa. Vseros konf. Tez. dokl. [Modern problems of remote sensing of the Earth from space. All-Russ. conf. Abstr.]. Moscow, 2005, p. 31.
11. Koviazin R.R. Dvumernoe vosstanovlenie fazy interferogramm [Two-dimensional reconstruction of interferograms phase]. Problemy kogerentnoi i nelineinoi optiki [Problems of coherent and nonlinear optics]. St. Petersburg, SPbGITMO Publ., 2000, pp. 267-275.
12. Filatov A.V. Metod obrabotki kompleksnykh radiolokatsionnykh interferogramm v usloviiakh vysokoi vremennoi dekorreliatsii. Kand.fiz.-mat. nauk diss. [Method of processing of complex radar interferograms at high temporal decorrelation. Cand. phys.-math. sci. diss.]. Barnaul, 2009.
13. El-taweel G.S. Enhanced model for generating 3-D images from RADAR interferometric satellite images. WSEAS Transactions on Environment and Development, 2007, vol. 3, no. 3, pp. 59-64.
14. Satalino G., Stramaglia S., Chiaradia M.T., Veneziani N. Weights determination for minimum cost flow InSAR phase unwrapping. Geoscience and Remote Sensing Symp. Proc. (IGARSS ’99). Hamburg, Germany, 1999, vol. 2, pp.1342-1344. DOI:10.1109/IGARSS.1999.774624
15. Ahuja R.K, Magnanti T.L., Orlin J.B. Network flows: Theory, algorithms, and applications. Englewood Cliffs, NJ, Prentice Hall, 1993.
