РАКЕТНО-КОСМИЧЕСКОЕ ПРИБОРОСТРОЕНИЕ И ИНФОРМАЦИОННЫЕ СИСТЕМЫ
2018, том 5, выпуск 1, c. 39–47
АЭРОКОСМИЧЕСКИЕ МЕТОДЫ ЗОНДИРОВАНИЯ ЗЕМЛИ
УДК 681.7.069.32+629.78
Обработка изображений на бортуКА ДЗЗ
методом инверсной фильтрации
В. Е.Квитка, аспирант МФТИ, kva-vasja@yandex.ru
Филиал АО РКЦ «Прогресс»–НПП «ОПТЭКС», г. Москва, Зеленоград, Российская Федерация
Московский физико-технический институт, г. Долгопрудный, Московская область, Российская Федерация
В. Д.Блинов, v_blinov@mail.ru
Филиал АО РКЦ «Прогресс»–НПП «ОПТЭКС», г. Москва, Зеленоград, Российская Федерация
Аннотация. Спутниковые снимки Земли высокого разрешения имеют размытие мелких деталей и малый контраст. Это обусловлено
низкими значениями функции передачи модуляции крупных оптических систем наблюдения Земли космического
базирования. В статье рассматриваются проблемы восстановления изображений методом инверсной фильтрации. Предлагаются
различные способы ее осуществления. Выявлена необходимость предварительного шумоподавления снимков перед их
восстановлением. Произведен количественный и качественный анализ работы наиболее известных алгоритмов шумоподавления
применительно к этой задаче. Разработаны и опробованы адаптивные алгоритмы предварительного шумоподавления.
Путем компьютерного моделирования изображений и сравнением критериев качества находится оптимальный для данной
задачи алгоритм предварительного шумоподавления снимков. Исследована возможность альтернативного решения проблемы
усиления шумов на снимке при выполнении инверсной фильтрации — применение регуляризации. Проведено компьютерное
моделирование различных способов регуляризации. Визуальный анализ результатов моделирования изображений и сравнение
количественных критериев их качества позволили определить оптимальные параметры функции регуляризации.
Ключевые слова: шумоподавление, инверсная фильтрация, функция передачи модуляции
Image Processing by Inverse Filtering
Onboard ERS Spacecraft
V.E. Kvitka, postgraduate student, kva-vasja@yandex.ru
Joint Stock Company “Space Rocket Centre «Progress»–Scientific and Production Enterprise «OPTEKS»”,
Zelenograd, Moscow, Russian Federation
Moscow Institute of Physics and Technology, Dolgoprudny, Moscow Region, Russian Federation
V.D. Blinov, v_blinov@mail.ru
Joint Stock Company “Space Rocket Centre «Progress»–Scientific and Production Enterprise «OPTEKS»”,
Zelenograd, Moscow, Russian Federation
Abstract. Satellite images of the Earth of high resolution have blurring of small details and low contrast. This is due to the
low values of the modulation transfer function of the large space-based optical Earth observation systems. The problems of image
reconstruction by the inverse filtering method are considered in the article. Various ways of its implementation are suggested.
The necessity of preliminary noise reduction of images before their restoration is revealed. A quantitative and qualitative analysis
of the performance of the most well known noise reduction algorithms has been performed. Adaptive algorithms of preliminary noise
reduction are developed and tested. By computer simulation of images and comparison of quality criteria, the optimal algorithm
of preliminary noise reduction of images is for this task is discovered. Regularization, an alternative solution to the problem of noise
amplification in an image when performing inverse filtering is studied. A computer simulation of various regularization methods
has been carried out. Visual analysis of the results of image modeling and comparison of quantitative criteria for their quality made
it possibleto determine theoptimal parameters of theregularization function.
Keywords: noise reduction, inverse filtering, modulation transfer function
Стр.1
40
В. Е.КВИТКА, В. Д. БЛИНОВ
Постановка задачи
Одно из свойств космических аппаратов дистанционного
зондирования Земли высокого разрешения
(КА ВР) — низкая величина функции
передачи модуляции (ФПМ) на высоких пространственных
частотах. Визуально это проявляется
в размывании изображений маленьких объектов
[1]. Исследуем возможность улучшения изображения
с помощью инверсной фильтрации, выполняемой
на борту космического аппарата.
Как известно, функция передачи модуляции
MTF(ν) связывает фурье-образ F(A) исходного
изображения A ифурье-образ F(B) изображения
B, получаемого оптико-электронной системой,
что показано уравнением (1). Здесь и далее речь
идет о двумерном дискретном преобразовании Фурье
(ДДПФ)
F(B)= F(A) · MTF(ν).
(1)
Естественным образом возникает предложение
восстановить изображение A,имеяизображение
B, полученное КА ВР при съемке поверхности
Земли. Проведем это действие по формуле (2):
A = F−1 F(B)
MTF(ν) .
(2)
Однако реализация данного метода, известного
как инверсная фильтрация, имеет большую
практическую трудность: при малых величинах
MTF(ν) шумы, всегда присутствующие на изображении
B, приведут к сильным искажениям восстановленного
изображения A [2]. Соответственно
возникает вопрос об исследовании границ применимости
данного метода восстановления изображений.
Ключевым моментом является необходимость
выполнения процедуры не на Земле, а при бортовой
обработке снимков. В противном случае будут
иметь место дополнительные искажения, связанные
с работой алгоритма сжатия видеоданных.
Моделирование снимков
Процедура моделирования снимков основана
на использовании аэроснимков в 12-битной градации
яркости. Так как рассматривается инверсная
фильтрация применительно к КА высокого разрешения
(около 0,5 м), то необходимо, чтобы исходные
снимки имели значительно меньшую проекцию пикселя.
В противном случае визуальный анализ результатов
инверсной фильтрации будет некорректным,
так как ФПМ аэрокамеры, существенно отличающаяся
от 1, внесет свои искажения. Поэтому
для моделирования используются аэрофотоснимки
с проекцией пикселя Laer не более 0,2 м.
Отметим его основные этапы:
1. Рассчитывается проекция пикселя аэрофотоснимка
Laer, после чего он перемасштабируется
на проекцию пикселя LКА (коэффициент масштабирования
k = Laer/LКА).Результатыэтой операции
описываются матрицей A.
2. Вычисляется функция передачи модуляции
(ФПМ) оптико-электронного тракта КА MTF(ν),
показанная на рис. 1.
3. Выполняется ДДПФ снимка A. Обозначим
результат как F(A).
4. Имитируем размытие мелких деталей изображения
A путем домножения спектра F(A) на
ФПМ оптико-электронного тракта MTF(ν).Полученный
результат F(B)= F(A) · MTF(ν) является
фурье-спектром изображения B,формируемого
в фокальной плоскости телескопа. Соответственно
для получения снимка B необходимо выполнить
обратное ДДПФ по формуле
B = F−1 F(B)) = F−1 F(A) · MTF(ν) .
5. Моделируем шумы снимка B сучетом фо(3)
тонного
шума на основе энергетического расчета
КА ВР для соответствующей высоты Солнца над
горизонтом. Тогда СКО шума имеет свое собственное
значение для каждого пикселя. Обозначим результат
как C.
Отметим, что здесь применена безразмерная
пространственная частота, нормированная на частоту
Найквиста (55,56 пар линий на мм).
На рис. 2–4 показаны фрагменты передискретизованного
аэрофотоснимка A, изображения на
фокальной плоскости объектива B и зашумленного
изображения C соответственно при возвышении
Солнца hs = 5◦. Все три снимка включают четыре
трехшпальных миры слева вверху (высокий и низкий
контрасты; периоды 2 и 4 пикселей).
РАКЕТНО-КОСМИЧЕСКОЕ ПРИБОРОСТРОЕНИЕ И ИНФОРМАЦИОННЫЕ СИСТЕМЫ т. 5 вып. 1 2018
Стр.2
ОБРАБОТКА ИЗОБРАЖЕНИЙ НА БОРТУ КА ДЗЗ МЕТОДОМ ИНВЕРСНОЙ ФИЛЬТРАЦИИ 41
Рис. 1. Функция передачи модуляции системы и ее аппроксимация
Рис. 2. Снимок A
Рис. 3. Снимок B
Исследование инверсной
фильтрации
Проведем инверсную фильтрацию без предварительного
шумоподавления для снимков C(hs),
зашумленность которых соответствует различным
возвышениям Солнца над горизонтом hs. Восстановленное
изображение Ainv описывается формулой
Ainv = F−1 F(C)
MTF(ν) .
(4)
Результат показан на рис. 5, где приведены
увеличенные фрагменты снимков.
Рис. 4. Зашумленный снимок C
Из рисунка видно, что удовлетворительный
результат инверсной фильтрации можно получить
только при возвышениях Солнца, превышающих
hs = 30◦, когда зашумленность снимков мала.
Тогда выполним сначала шумоподавление, получив
снимок Cfil, и только после этого проведем инверсную
фильтрацию над ним. Результаты показаны
на рис. 6.
Из рисунка видно, что лучший результат обеспечивается
при инверсной фильтрации снимка, прошедшего
шумоподавление псевдомедианным фильтром,
описанным в [3]. В этом методе сигнал пикселя
заменяется на med (y1, y2, y3),где y1, y2,
y3 — значения медиан сигнала в строках окна
РАКЕТНО-КОСМИЧЕСКОЕ ПРИБОРОСТРОЕНИЕ И ИНФОРМАЦИОННЫЕ СИСТЕМЫ т. 5 вып. 1 2018
Стр.3
42
В. Е.КВИТКА, В. Д. БЛИНОВ
Рис. 5. Инверсная фильтрация без предварительного шумоподавления
Рис. 6. Результаты шумоподавления и последующей инверсной фильтрации
3 Ч 3 пикселя. Восстановление после медианной
фильтрации шума породило искажения трехшпальной
миры, а среднегеометрический фильтр настолько
сильно размывает детали, что инверсная фильтрация
не способна восстановить резкость краев.
Предварительная обработка
адаптивным фильтром
Исследуем инверсную фильтрацию, выполненную
над снимком, прошедшим шумоподавление
адаптивным алгоритмом. Идея такого фильтра описана
в [2, с. 355]. Он должен учитывать зависимость
СКО шума от сигнала в пикселе, а при наличии
сильных перепадов яркости не сглаживать
их, так как это не шумы, а контуры объектов.
Предлагается следующий алгоритм: в качестве
признака наличия контуров в окне 3 Ч 3пикселя
E примем условие, что СКО2(E)/noise2(i, j)9.
То есть фактическое СКО шума в окне 3 Ч 3пикселя
в 3 раза больше, чем шум в пикселе C(i, j).
Блок-схема алгоритма показана на рис. 7.
При этом СКО шума определяется бортовым
устройством по формуле
noise(i, j)= 4095
Pзар · C(i, j) · Pзар
4095.
(5)
Зарядовая емкость фотоприемника Pзар известна
заранее. Однако вычисление квадрата величины
СКО сигнала в окне 3 Ч 3пикселя требуетоколо
Выявление контуров с помощью
градиента
Используем для шумоподавления градиентный
алгоритм, учитывающий направление контуроввскользящем
окне3×3 пикселя. Как известно,
РАКЕТНО-КОСМИЧЕСКОЕ ПРИБОРОСТРОЕНИЕ И ИНФОРМАЦИОННЫЕ СИСТЕМЫ т. 5 вып. 1 2018
20 операций, что весьма много для задачи бортовой
обработки.
Стр.4
ОБРАБОТКА ИЗОБРАЖЕНИЙ НА БОРТУ КА ДЗЗ МЕТОДОМ ИНВЕРСНОЙ ФИЛЬТРАЦИИ 43
Рис. 7. Блок-схема алгоритма
градиент функции двух переменных U(x, y) является
вектором, задающим направление наибольшего
изменения функции. При наличии в окне контура
вектор-градиент будет перпендикулярен резкому
краю. Тогда усредним яркости пикселей именно
по направлению контура, чтобы не допустить
его размывания. Математически это описывается
формулой
Cfil(i, j)= a · C(i, j)+
+ b · sin2 θ · C(i, j − 1)+C(i, j + 1)
2
2
+
+cos2 θ · C(i − 1, j)+C(i + 1, j)
. (6)
При этом выполняется условие a+b = 1, необходимое
для неизменности средней яркости изображения
однородной поверхности. Угол θ, являющийся
углом между вектором-градиентом и осью X,
определяется по формуле
tg θ =
∂U
∂y
∂U
∂x
= C(i + 1, j) −C(i, j)
C(i, j + 1) −C(i, j) .
Тогда sin2 θ = tg2 θ
1+tg2 θ,а cos2 θ = 1
1+tg2 θ
.
Эти равенства позволяют избежать трудоемкого
вычисления синуса и косинуса методом разложения
в ряд Тейлора.
Отметим, что угловые пиксели окна 3 Ч 3
не участвуют в вычислениях. Для улучшения передачи
негоризонтальных и невертикальных контуров
введем в алгоритм ветвь, учитывающую величины
угловых пикселей окна 3 Ч 3. Блок-схема
(7)
градиентного алгоритма шумоподавления показана
на рис. 8.
При рассмотрении диагонального контура
30◦ < |θ| < 60◦ одна ветвь работает, когда контур
можно считать диагональным и при этом в верхней
части окна 3×3 контур находится вверху, а в правой
части окна — внизу. Аналогично другая ветвь
работает, если в верхней части окна 3 Ч 3диагональный
контур находится внизу. Результаты применения
алгоритма для значений параметров a =
= 0,5 и b = 0,5 показаны на рис. 9–11.
Как видно из рис. 9–11, при высоте Солнца,
превышающей 10◦, обеспечивается вполне приемлемый
результат инверсной фильтрации.
Шумоподавление на основе
лапласиана
Выявлять контуры и неоднородности в окне
3 Ч 3 можно с помощью оператора Лапласа (лаплаLap(i,
j)= ∂2U
∂x2 + ∂2U
∂y2 =
= (C(i, j − 1)+C(i, j + 1) − 2 · C(i, j))+
+(C(i − 1, j)+C(i + 1, j) − 2 · C(i, j)) . (8)
Очевидно, что для белой точки на черном фоне
получится, что Lap(i, j)= 2 ·|0 + 0 − 4095| = 8160.
Для абсолютно резкой горизонтальной белой линии
на черном фоне имеем, что Lap(i, j)= |0 + 0 −
− 4095| = 4095.
РАКЕТНО-КОСМИЧЕСКОЕ ПРИБОРОСТРОЕНИЕ И ИНФОРМАЦИОННЫЕ СИСТЕМЫ т. 5 вып. 1 2018
сиана). Определим его в уравнении
Стр.5
44
В. Е.КВИТКА, В. Д. БЛИНОВ
Рис. 8. Блок-схема градиентного алгоритма шумоподавления
Рис. 9. a = 0,5, hs = 2◦
Рис. 10. a = 0,5, hs = 5◦
Рассмотрим алгоритм, основанный на трех пороговых
значениях лапласиана в текущей точке
снимка C(i, j). Блок-схема показана на рис. 12.
Здесь Sij —окно3 Ч 3 пикселя, а Sij —
проколотая окрестность точки (i, j),тоестьокно
3 Ч 3 без центрального пикселя. Соответственно
Рис. 11. a = 0,5, hs = 10◦
Применим описанный алгоритм для предварительного
шумоподавления, проводимого перед инверсной
фильтрацией. Результат получения снимков
Ainv для разных возвышений Солнца показан
на рис. 13–14.
Как видно из рис. 13–14, результат восстааSij
обозначает среднее значение по 9 пикселям,
Sij — среднее значение сигнала 8 пикселей.
Таким образом, мы отказываемся от изменения
сигнала в точках с очень сильными локальными
неоднородностями (Lap(i, j) 1000), чтобы
избежать размытия мелких деталей. Наоборот, для
равноярких участков работает ветвь, усредняющая
с одинаковыми весами все 9 пикселей окрестности.
новления
для данного алгоритма предварительного
шумоподавления хуже, чем при использовании градиентного
метода, так как имеют место одиночные
темные или светлые точки на однородных участках.
Регуляризация решения
Кроме применения адаптивных алгоритмов
шумоподавления, проблему деления фурье-образа
РАКЕТНО-КОСМИЧЕСКОЕ ПРИБОРОСТРОЕНИЕ И ИНФОРМАЦИОННЫЕ СИСТЕМЫ т. 5 вып. 1 2018
Стр.6
ОБРАБОТКА ИЗОБРАЖЕНИЙ НА БОРТУ КА ДЗЗ МЕТОДОМ ИНВЕРСНОЙ ФИЛЬТРАЦИИ 45
Рис. 12. Блок-схема лапласианного алгоритма шумоподавления
Рис. 13. Ainv, hs = 2◦
Рис. 14. Ainv, hs = 10◦
на малую величину ФПМ можно решить регуляризацией
решения, описанной в [4]. Идея регуляризации
описывается формулой
Ainv = F−1
MTF(ν)+Reg(ν) ,(9)
F(C)
где Reg(ν) — функция регуляризации.
Из анализа формулы (9) можно сформулировать
два требования к функции регуляризации:
ных • Reg(0)= 0, так как на малых пространствена
к 1, чтобы не допускать большого искажения
частотах искажения минимальны и регуляризация
не нужна;
к • Reg(1) должна быть ненулевая, но и не близвосстановленного
изображения.
Использование регуляризации вместо предварительного
шумоподавления представляет интерес,
так как требует меньшего числа операций с массивом
данных. Проведем анализ различных вариантов
функции регуляризации.
Квадратичная зависимость
Reg(ν)= α · ν2,где α — экспериментально
подбираемый безразмерный коэффициент. Попытаемся
экспериментально оптимизировать величину
α. На рис. 15 показано моделирование инверсной
фильтрации для разных значений α без предварительного
шумоподавления.
При анализе рис. 15, а–в видно, что с ростом α
возрастает соотношение сигнал/шум, но при этом
падает передача контраста. Характер данной зависимости
показан на рис. 16. Далее будем говорить об
отношении сигнал/шум, приведенном для сигнала,
составляющего половину от максимума, что соответствует
альбедо 0,5. Обозначим данную величину
SNR50.
РАКЕТНО-КОСМИЧЕСКОЕ ПРИБОРОСТРОЕНИЕ И ИНФОРМАЦИОННЫЕ СИСТЕМЫ т. 5 вып. 1 2018
Стр.7
46
В. Е.КВИТКА, В. Д. БЛИНОВ
Рис. 15. Моделирование инверсной фильтрации без предварительного шумоподавления
Рис. 16. Зависимости коэффициента передачи контраста (КПК) и изменения соотношения сигнал/шум
(SNR после фильтрации/SNR исходное) от коэффициента α
Принимая, что КПК должен быть близок к 1,
как и изменение SNR, определим, что оптимальное
значение α находится в пределах от 0,6 до 1.
Кубическая зависимость
Рассмотрим другую регуляризующую функцию
— Reg(ν)= β · ν3,где β —эксперименполучаем,
что оптимальное значение β находится
в пределах от 0,8 до 1,2.
Заключение
Для математического описания качества инверсной
фильтрации сведем в табл. основные хатально
подбираемый безразмерный коэффициент.
Тогда влияние изменения знаменателя дроби (3)
на малыхчастотахбудет меньше,чем уквадратичной
регуляризации. На рис. 15, а–в показано моделирование
инверсной фильтрации без шумоподавления.
Аналогично квадратичной регуляризации
сростом β возрастает соотношение сигнал/шум,
но при этом падает КПК. Характер данной зависимости
показан на рис. 17.
Принимая, что коэффициент передачи контраста
должен быть близок к 1, как и изменение SNR,
рактеристики ее вариантов. Здесь ФПК (0,5νN)—
отношение контрастов миры с пространственной
частотой 50% частоты Найквиста для преобразованного
и исходного снимков. Отрицательное
значение ФПК соответствует инверсии шпал
миры: черная шпала становится белой и наоборот.
Аналогичное явление можно наблюдать на рис. 6
при псевдомедианной фильтрации.
Проанализируем данные таблицы и результаты
работы предварительных алгоритмов шумоподавления.
Оптимальными решениями для подготовки
РАКЕТНО-КОСМИЧЕСКОЕ ПРИБОРОСТРОЕНИЕ И ИНФОРМАЦИОННЫЕ СИСТЕМЫ т. 5 вып. 1 2018
Стр.8
ОБРАБОТКА ИЗОБРАЖЕНИЙ НА БОРТУ КА ДЗЗ МЕТОДОМ ИНВЕРСНОЙ ФИЛЬТРАЦИИ 47
Рис. 17. Зависимости КПК и изменения соотношения сигнал/шум (SNR после фильтрации/SNR исходное)
от коэффициента β
Та б лиц а. Основные характеристики вариантов инверсной фильтрации
Вариант предварительного
шумоподавления
Без шумоподавления
Псевдомедианой
Градиентный алгоритм
Лапласианный алгоритм
Регуляризация α · ν2 (α = 0,8)
Регуляризация β · ν3 (β = 1)
SNR50
31
31
31
31
31
изображений представляются использование градиентного
метода шумоподавления или же кубическая
регуляризация. Кроме того, можно применять
совместное использование этих операций, когда регуляризация
выполняется для фурье-образа снимка,
прошедшего шумоподавление.
Список литературы
1. Бакланов А. И. Системы наблюдения и мониторинга:
учебное пособие. М.: БИНОМ. Лаборатория знаний,
2009. 234 с.
5,6
15
9
21
31 · 0,7 = 22
31 · 0,75 = 23
ФПК
на исходном снимке на обработанном снимке 0,5νN νN
31
3
7
8,26
2,4 −5
3
3,7
1
1
8,6
1
1
2. Гонсалес Р., Вудс Р. Цифровая обработка изображений.
М.: Техносфера, 2005. 1072 с.
3. Акиншин Н. С., Есиков О. В.Псевдомедианная фильтрация
для обработки оптических изображений //
Материалы XII научно-технической конференции
«Системы наблюдения, мониторинга и дистанционного
зондирования Земли». Калуга: Манускрипт, 2015.
436 с.
4. Кашкин В. Б. Цифровая обработка аэрокосмических
изображений. Версия 1.0 [Электронный ресурс]: конспект
лекций / В. Б. Кашкин, А.И. Сухинин. Красноярск:
ИПК СФУ, 2008.
РАКЕТНО-КОСМИЧЕСКОЕ ПРИБОРОСТРОЕНИЕ И ИНФОРМАЦИОННЫЕ СИСТЕМЫ т. 5 вып. 1 2018
Стр.9