CC BY
13
5. Weizsäcker С. С., Ferner J. An integrated simulation model for European electricity and natural gas supply // Electrical Engineering. 2001. 83. N 5. P. 265-270.
6. Гантмахер Ф.Р. Теория матриц. M.: Наука, 1966.
7. Ланкастер П. Теория матриц. М.: Наука, 1973.
8. Каре лова О. Л., Бань ко М. А. Применение марковских цепей при прогнозировании демографической ситуации в мире // Математическое моделирование. 2006. 18. № 2. С. 43-50.
9. Выдрин A.C., Михалёв A.B. Стохастические матрицы и анализ защищенности автоматизированных систем // Фундаментальная и прикладная математика. 2007. 13. № 1. С. 61-99.
10. James С. Fu, Liqun Wang, Wendy Lou W.Y. On exact and large deviation approximation for the distribution of the longest run in a sequence of two-state Markov dependent trials //J. Appl. Probability. 2003. 40. N 2. P. 346-360.
11. Seneta E. Nonnegative Matrices and Markov Chains. N.Y.: Springer, 1981.
Поступила в редакцию 05.10.09
УДК 519.612:632.4 В.А. Крылов1
АППРОКСИМАЦИЯ РАСПРЕДЕЛЕНИЙ АМПЛИТУД ИЗОБРАЖЕНИЙ РАДАРА С СИНТЕЗИРОВАННОЙ АПЕРТУРОЙ МЕТОДОМ КОНЕЧНЫХ СМЕСЕЙ
В работе предлагается метод для решения проблемы аппроксимации распределения амплитуд для изображений, полученных радаром с синтезированной апертурой (РСА), с использованием метода конечных смесей распределений. Данный метод включает в себя стохастический Expectation Maximization (ЕМ) алгоритм и метод логарифмических кумулянт для оценки параметров компонент смеси. Распределения для компонент смеси выбираются из специального словаря, содержащего характерные для РСА распределения. Предложенный метод является полностью автоматическим. Эксперименты с реальными РСА-изображениями высокого разрешения показывают результаты высокой точности как с точки зрения визуального анализа, так и с точки зрения количественных характеристик (коэффициент корреляции, расстояние Колмогорова-Смирнова).
Ключевые слова: радар с синтезированной апертурой, распределение амплитуд, метод конечных смесей, стохастический ЕМ-алгоритм, метод логарифмических кумулянт.
1. Введение. В контексте задач дистанционной регистрации данных одна из важнейших задач связана с аппроксимацией распределений амплитуд изображений, полученных со спутников и радаров, в частности при помощи радаров с синтезированной апертурой (РСА) [1]. Приближенные распределения амплитуд используются в ряде дальнейших прикладных задач, таких, как удаление шумов, классификация и распознавание объектов. В литературе предложен ряд распределений, хорошо аппроксимирующих наблюдаемые распределения амплитуд РСА, эти распределения можно условно разбить на два класса: эмпирические и теоретически-обоснованные. К первому классу относятся: ло-гнормальное распределение, распределение Вейбулла и обобщенное гамма-распределение [1]; эти распределения, не имея теоретических обоснований применимости к РСА, показали хорошие результаты в экспериментах. В свою очередь, для теоретически-обоснованных распределений, таких, как распределение Накагами,
К1/2
-распределение (К-распределение описывает распределение интенсивностей, а соответствующее распределение К1/2 было предложено для амплитуд [1]) и обобщенное распределение Гаусса-Рэлея [1, 2], были предложены доказательства (с дополнительными предположениями) соответствия этих моделей определенным видам поверхностей (вода, лес, поле) [1]. Однако слабость всех приведенных выше моделей заключается в невозможности описания распределений неоднородных
1 Факультет ВМК МГУ, асп., e-mail: vkrylovQcs.msu.ru
сцен, где на гистограмме одновременно присутствуют характерные особенности нескольких различных видов поверхностей. Для решения этой проблемы в работе предлагается метод конечных смесей распределений, где компоненты смеси берутся из приведенных выше распределений.
2. Метод конечных смесей и словарь. Имеем амплитудное РСА-изображение в виде множества X = {г!,г2,... ,гдг} независимых, одинаково распределенных наблюдений (н.о.р.н.) из конечной смеси независимых распределений:
к
= г^ о, (1)
1=1
к
где Рг(г|$г) — параметризованные плотности с 0, (г (-),• С . {1\.....Рк} — веса смеси: V I) I и
г=1
О ^ Рг ^ 1, г = 1,..., .К". Задача состоит в оценке вектора в = (01,..., 9к; Ръ • • •, Рк)- Подход н.о.р.н., игнорируя контекстную зависимость наблюдений, позволяет применить метод конечных смесей.
Каждая из плотностей в (1) выбирается из словаря V = {/1,..., /а/}, состоящего из М = 6
(см. табл. 1) характерных для РСА плотностей распределений /¿(ж|аг). В табл. 1 также приведены уравнения метода логарифмических кумулянт, описанного в п. 4. Подробное описание свойств используемых распределений можно найти в [3].
Табл и ца 1
Распределения, составляющие словарь V*
Распределение Плотность MJIK-уравнения
Логнормальное Kl = rn, К2 = а2
Вейбулла h(r) = JLrV-lexp Ki = In fj, + Ф(1)?7"1, К2 = Ф(1,1)?Г2
Обобщенное гамма (ОГ) Л(г) = ^fo (¿Г"1 «р {-(£)"} Kl = + In <7, К2 = Ф(1 ,K)/v2, КЗ = Ф(2,
Накагами 2кг = In fj, + Ф(L) - In L, АК2 = 9(1, L)
Ki/2 h(r) = Г(Ь)Г(М) х Хгь+м-1 Км-L [2r (LMm)1-72] 2KI = Ф(L) + Ф(М) - In(LMfi), 4К2 = Ф(1, L) + Ф(1,М), 8к3 = Ф(2, L) + Ф(2, М)
Обобщенное Гаусса-Рэлея (ОГР) /б(г) = -^Фщ /2«p[-(7r)1/A«W]d0, s(0) = |cos6>|1/A + |sin6>|1/A KI = АФ(2А) - 1п7--AGi(A)[G0(A)]-1, кг = -A2[Gi(A)]2[G0(A)]-2 + +А2Ф(1, 2А) + A2G2(A)[G0(A)]-1
* Г(-) — гамма функция [7], Ка(-) — модифицированная функция Бесселя второго рода порядка а [7], Ф(-) — дигамма функция, •) — полигамма функция порядка V [3], Ср(-) — интегральные функции, определяющие ОГР [2].
3. Стохастический ЕМ-алгоритм. Алгоритм Expectation Maximization (ЕМ) является одним из классических методов локальной оптимизации, в частности, он активно применяется в задачах оценки смесей [4]. ЕМ-алгоритм решает задачу с неполными данными, где по неполному вектору наблюдений у € Y требуется восстановить полный вектор х = (у, z) € X. В терминологии смесей для каждого наблюдения у^ требуется восстановить номер компоненты Zj, которой принадлежит данное наблюдение. ЕМ-алгоритм локальной оптимизации крайне чувствителен к инициализации. Одним из способов ухода от проблемы инициализации является использование одной из модификаций ЕМ-алгоритма — стохастического ЕМ (SEM) [5]. Начиная с произвольной начальной конфигурации, в SEM-алгоритме итеративно повторяются следующие шаги:
- Е-шаг: вычисление условной плотности рх\у(-\у, в1) с применением текущей оценки вг;
- S-шаг: построение реализации полного вектора ж* € X по плотности, полученной на Е-шаге;
- М-шаг: получение оценки максимального правдоподобия (ОМП) 9t+l € © по реализации ж*, построенной на S-шаге.
Последовательность SEM-оценок является дискретным случайным процессом, который не сходится ни поточечно, ни почти наверное, однако является однородной марковской цепью [5]. В случае конечных смесей SEM-процесс также является эргодической цепью, что обеспечивает сходимость к стационарному распределению, "сконцентрированному" в окрестности точки максимума функции максимального правдоподобия (иными словами, оценка максимального правдоподобия параметров смеси является асимптотически эквивалентной математическому ожиданию вг относительно стационарного распределения SEM-процесса).
Благодаря S-шагу алгоритм "выходит" из неустойчивых локальных максимумов, оставаясь, тем не менее, алгоритмом локальной оптимизации. Дополнительное достоинство SEM состоит в возможности применять его и для оценки количества компонент в смеси [5]. Идея состоит в инициализации алгоритма оценкой сверху для количества компонент, позволяя по ходу работы алгоритма исключать компоненты, вес которых становится слишком малым. Эксперименты в [5] показывают хорошие результаты при использовании такой оценки размера смеси. Для выбора начальной точки х° € X алгоритма предлагается применение поиска мод эмпирического распределения с последующей инициализацией нескольких компонент в окрестности каждой моды гистограммы. Такой подход можно обосновать использованием словаря исключительно из одномодальных распределений [1]. Такая инициализация также удовлетворяет условию оцененного сверху начального количества компонент.
4. Метод логарифмических кумулянт. В классической схеме SEM на М-шаге предлагается использование ОМП, однако, для К1//2 из V поиск этой оценки вычислительно слишком затруднителен [1]. Поэтому для оценки параметров на М-шаге предлагается использование метода логарифмических кумулянт (МЛК) [6]. Этот метод был разработан по аналогии с классическим методом моментов (ММ), где оценка параметров находится из системы уравнений, связывающих аналитические выражения для моментов и их эмпирическую оценку. В [6] ММ адаптируется для с.в. с плотностью на [О,+оо) (как распределения амплитуд РСА), используя по аналогии с преобразованием Лапласа в производящей функции моментов интегральное преобразование Меллина [7]:
+ ОС
</>u(s) = J pu{u)us~l du, se С. о
Определяются кумулянты и-го порядка как к„ = [1ik/>„(s)]^(1), v = 1, 2,..., где ^ означает производную v-то порядка. В случае сходимости преобразования имеют место следующие МЛК-уравнения [6]:
Ki = Е{Ыи}, «2 = Var{lnti}, кз = Е{(Ыи — Ki)3}.
Как видно из табл. 1, системы МЛК-уравнений в некоторых случаях не имеют аналитического решения относительно искомых параметров, однако благодаря строгой монотонности фигурирующих в уравнениях функций, численное обращение этих систем не представляет сложности. Стоит также отметить, что МЛК-оценки имеют в некоторых случаях хорошие теоретические свойства, такие, как, например, состоятельность в случае ОГР [2].
5. Общая структура алгоритма. После инициализации, описанной в п. 3, итерационно повторяются следующие шаги:
- Е-шаг: вычисление для каждого уровня в гистограмме г и компоненты г апостериорной вероятности в соответствии с текущей оценкой, 0.....X - \. / I.....К:
тЦг) = кРЫ(г) ; £
- 8-шаг: реализация меток я1 (г) каждого элемента гистограммы г в соответствии с мультиномиальными распределениями (т/(г): г = 1,..., К), 0.....7. — \\
- К-шаг: удаление таких компонент й, й = 1 что Р\ менее некоторого порогового значения; обновление числа компонент К:
- МЛК-шаг: для каждой компоненты г, вычисление веса компоненты и кумулянт:
£ Кг) £ Цг) Ыг £ Цг)(Ыг - к\г)ь
1 =
Z-l
Е ВД
2 = 0
I _ zeQit
£ М*)
zeQit
I геЯи
"5 КЫ —
£ М*)
zeQit
где Ь = 2,3, 1г(г) — гистограмма, (¿ц = : (г) = сг^} — уровни гистограммы, входящие в г-ю компоненту; затем решение системы МЛК-уравнений (см. табл. 1);
- шаг выбора модели: для каждой компоненты г, вычисление логарифма функции правдоподобия
Л-('К-)
ЬЬ = Кг)Ъь1М<*\з)
и выбор в качестве плотности имеющей наибольшее значение г = 1 ,
3 = 1 ,...,м.
6. Эксперименты. Для демонстрации выбраны несколько изображений, полученных немецким РСА ТеггаЭАИ-Х. Изображения РСА1 и РСА2 (рисунок) содержат неоднородные сцены и соответствующие им гистограммы многомодальны, что указывает на низкую точность при аппроксимации одномодальным распределением. Эти изображения получены с разными параметрами РСА, однако имеют высокое разрешение: 6 метров на пиксель для РСА1 и около 2.5 для РСА2 и РСАЗ.
Изображения РСА1 (а) и РСА2 (в) и соответствующие им гистограммы: исходная, смеси, отдельных компонент смеси и лучшего одномодельного приближения из Т> для РСА1 (б) и для РСА2 (г)
Количественная оценка аппроксимации приводится в табл. 2. Для сцен РСА1 и РСА2 повышение точности существенно при адекватном количестве компонент К = 3. Сцена на изображении РСАЗ практически однородная и, как видно из табл. 2, улучшение аппроксимации незначительно. Эта ситуация характерна для построенной модели: ее преимущество проявляется при аппроксимации многомодальных гистограмм, соответствующих неоднородным сценам.
Табл и ца 2
Количественная оценка аппроксимации смесью и лучшей моделью из V*
Изображение Смешанная модель Лучшая модель из D
КС К Смесь Модель КС
РСА1 0.008 3 (/3,/4,/х) Накагами 0.037
РСА2 0.011 3 (h,h,fi) ОГР 0.029
РСАЗ 0.007 3 (/4) /з) /4) Вейбулл 0.020
* КС — расстояние Колмогорова-Смирнова между приближением и исходной гистограммой, К — оценка количества компонент.
Отдельно отметим, что благодаря способности моделировать распределение неоднородных классов, разработанный метод представляется перспективным для применения к задаче классификации поверхностей на PC А. В стандартной постановке этой задачи проводится спецификация интересующих классов поверхностей, каждый из которых зачастую является смесью нескольких подклассов поверхностей, что может быть аккуратно описано в рамках предложенной модели.
7. Заключение. В работе представлен автоматический алгоритм для моделирования распределения амплитуд PC А-изображений на основе использования метода конечных смесей. Для оценки модели применяются методы стохастического ЕМ и логарифмических кумулянт. Эксперименты на реальных изображениях высокого разрешения показывают высокую точность аппроксимаций, построенных разработанным методом. Метод представляется перспективным для применения к задаче классификации поверхностей на РСА изображениях.
Автор выражает благодарность своему научному руководителю доценту Матвееву В. Ф., а также профессорам Мозеру Г., Серпико С. и Зерубии Д. за помощь при подготовке этой работы.
Использованные в работе РСА-изображения системы TerraSAR-X распространяются на свободной основе (http://www.infoterra.de/, ©Infoterra).
СПИСОК ЛИТЕРАТУРЫ
1. Oliver С., Quegan S. Understanding Synthetic Aperture Radar Images. Norwood: Artech House, 1998.
2. Moser G., Zerubia J., Serpico S.B. SAR amplitude probability density function estimation based on a generalized gaussian model // IEEE Transactions on Geoscience and Remote Sensing. 2006. N 15. P. 1429-1442.
3. Krylov V., Moser G., Serpico S.B. et al. Modeling the statistics of high resolution SAR images. Research Report. N 6722. Sophia Antipolis: INRIA, 2008.
4. Red ner R. A., Walker H. F. Mixture densities, maximum likelihood, and the EM algorithm // SIAM Review. 1984. 26. N 2. P. 195-239.
5. Celeux G., Chauveau D.,Diebolt J. On stochastic versions of the EM algorithm. Research Report. N 2514. Grenoble: INRIA, 1995.
6. Nicolas J.-M., Maruani A. Lower-order statistics: a new approach for probability density functions defined on R+ // Proc. European Signal Processing Conf. EUSIPCO 2000. Tampere: EURSICO, 2000. P. 805-808.
7. Sneddon I. The Use of Integral Transforms. N.Y.: McGraw-Hill, 1972.
Поступила в редакцию 25.11.09
