Математическое моделирование процесса радиометрической коррекции снимков дистанционного зондирования Земли Текст научной статьи по специальности «Математика»

УДК 51.74, 528.854, 528.852

Математическое моделирование процесса радиометрической коррекции снимков дистанционного зондирования Земли

© В.И. Майорова, А.М. Банников, К.И. Зайцев МГТУ им. Н.Э. Баумана, Москва, 105005, Россия

Представлена математическая модель процесса радиометрической коррекции изображения. Рассмотрены методы определения искажающей функции системы. Приведены и программно реализованы следующие алгоритмы восстановления изображений: инверсная фильтрация, фильтрация методом минимизации среднеквадратического отклонения (фильтрация Винера), фильтрация методом минимизации сглаживающего функционала со связью (регуляризация Тихонова). Проведено сравнение перечисленных алгоритмов, выявлены их преимущества и недостатки применительно к восстановлению космических снимков.

Ключевые слова: искажающая функция, оптическая передаточная функция, инверсная фильтрация, фильтрация методом минимизации среднеквадратического отклонения (фильтрация Винера), фильтрация методом минимизации сглаживающего функционала со связью (регуляризация Тихонова), функционал невязки.

Введение. В настоящее время актуальна проблема разработки новых космических средств дистанционного зондирования Земли (ДЗЗ) для их применения во многих отраслях народного хозяйства. Неотъемлемым элементом средств ДЗЗ является датчик — электронно-оптическое устройство, преобразующее падающее на него отраженное солнечное излучение или собственное излучение Земли в цифровой снимок с пространственным распределением излучения. Однако вследствие несовершенства конструкции датчика, оптических эффектов, электронных шумов, бега подстилающей поверхности и других внутренних и внешних факторов ухудшается качество принимаемого сигнала, т. е. происходит искажение информации о земной поверхности. Математическое моделирование этого искажения помогает уяснить суть процесса появления искажений, что необходимо для понимания самой методики и работы алгоритмов радиометрической коррекции космических снимков.

В статье рассмотрено математическое моделирование радиометрической коррекции снимков ДЗЗ с использованием искажающей функции. Эта функция устанавливает взаимосвязь между идеальным (неискаженным) изображением и реальным (искаженным) изображением на выходе электронно-оптической системы.

Для моделирования разработаны программы в среде Ма1ЪаЬ. Проведены эксперименты по восстановлению смазанных, размытых и зашум-ленных реальных снимков ДЗЗ, выполненных с космического аппарата «Landsat-5» с применением таких алгоритмов восстановления цифровых изображений, как инверсная фильтрация, фильтрация методом минимизации среднеквадратического отклонения (фильтрация Винера), фильтрация методом минимизации сглаживающего функционала со связью (регуляризация Тихонова). Получены результаты восстановления снимков. Выявлены преимущества и недостатки каждого из перечисленных алгоритмов применительно к процессу восстановления снимков ДЗЗ.

Математическое моделирование процесса радиометрической коррекции изображения. Процесс искажения изображения в общем

виде можно описать следующей формулой:

g(x, у) = H [ / & у)] + у\

где g(x, у) — искаженное изображение (физический сигнал на выходе оптической системы); /(х, у) — идеальное (неискаженное) аналоговое изображение; Н — искажающий оператор оптической системы, п (х, у) — аддитивный шум.

Отметим, что функции g(x, у), /(х, у), п(х, у) являются непрерывными функциями от непрерывных пространственных координат, поскольку преобразование к дискретным величинам произойдет только после прохождения процедуры аналого-цифрового преобразования сигнала.

Имея функцию g(x, у), обладая информацией об искажающем операторе Н оптической системы и зная основные характеристики аддитивного шума п^, у), можно построить некоторое оптимальное приближение /(х, у). При этом чем больше имеется информации об операторе Н и аддитивном шуме п^, у), тем точнее может быть построено

приближение / (х, у). Схема, моделирующая процесс искажения и коррекцию изображения, представлена на рис. 1.

Если выполняются условия

• линейности Ща/^, у) + Ъ^, у)) = аН( у^, у)) + ЪН( у)) и

• трансляционной инвариантности Н(у ^ - а, у - Р)) = g(x - а, у - Р), можно показать математически, что искаженное изображение предста-вимо в пространственной области в следующем виде [1]:

g(x, у) = у) * /у) + п(^ уХ (1)

где И^, у) — искажающая функция (пространственное представление искажающего оператора Н ). В различных источниках функцию И^, у) называют также функцией пространственного отклика, ядром иска-

Рис. 1. Моделирование процесса искажения и коррекции изображения

жающего оператора, функцией kernel, функцией рассеяния точки (PSF, от англ. Point Spread Function) [2-3].

Операция применения искажающей функции к другой функции (в данном случае к изображению f (x, y)) называется сверткой (от англ. convolution). В формуле (1) эта операция обозначена символом *. Смысл операции свертки состоит в том, что в формировании значения точки A с координатами (x0, y0), принадлежащей искаженному изображению g, принимают участие точки из некоторой окрестности 5 точки B с координатами (x0, y0), принадлежащие исходному изображению f причем размер окрестности 5 и вклад точек из окрестности B в формирование значения точки A целиком определяется искажающей функцией.

Для пояснения вышесказанного рассмотрим пример. Предположим, что ведется работа с растровым (пиксельным) изображением, имеющим разрешение M х N. В таком случае переменные x, y будут определять порядковый номер пикселя в изображении, следовательно, x е 0, 1, ...,M; y е 0, 1, ..., N. Пусть окрестность 5 имеет размер m х n. Тогда математически процесс применения операции свертки к изображению записывается следующим образом:

ab

g(х, y) = h(x, y )* f (x, y ) = — £ £ h(i, j) f (x + i, у + j), (2)

m -1 ; ь = n -1

2 2

На практике часто а = Ь. В этом случае число а называется размерностью искажающей функции. Размерность искажающей функции, как правило, меньше размерности самого изображения.

На рис. 2 представлен пример искажающей функции для моделирования оптического размытия изображения, а на рис. 3 — пример искажающей функции для моделирования смазывания изображения (размерность искажающей функции в обоих примерах равна 30).

В соответствии с теоремой о свертке, в пространственной области свертка функций эквивалентна умножению в частотной области Фурье-

Рис. 2. Искажающая функция размерности 30 для моделирования размытия по закону Гаусса

Рис. 3. Искажающая функция для моделирования смазывания изображения за счет бега подстилающей поверхности

преобразований этих функций, поэтому приведенное выше уравнение модели искажения в пространственной области (2) можно записать в эквивалентном представлении в частотной области:

G(u, v) = H(u, v)F(u, v) + N(u, v),

(3)

где заглавными буквами даны соответствующие Фурье-преобразования функций, обозначенных строчными буквами в уравнении свертки в пространственной области (2).

Функцию H(u, v) часто называют оптической передаточной функцией (OTF, от англ. Optical Transfer Function). Этот термин заимствован

из анализа Фурье оптических систем. Очевидно, что функции h(u, v) и H(u, v) переходят одна в другую под действием прямого и обратного Фурье-преобразования.

Вследствие того, что оператор линейного трансляционно-инвари-антного искажения Н можно смоделировать в виде свертки (конволю-ции), иногда этот процесс искажения называют конволюцией изображения с PSF или OTF. По аналогии обратный процесс восстановления называется деконволюцией.

Моделирование искажающей функции. Получение точной оценки искажающей функции оптико-электронной системы — один из важных шагов на пути к восстановлению искаженного изображения. Существует три основных способа оценки искажающей функции оптико-электронной системы [4]:

• визуальный анализ;

• экспериментальный анализ;

• математическое моделирование.

Рассмотрим каждый из этих способов.

Визуальный анализ. Предположим, что имеется искаженное изображение, но информация об искажающей функции Н отсутствует. Один из способов оценить эту функцию состоит в выделении информации непосредственно из изображения. Например, если изображение является размытым, можно рассмотреть его небольшой фрагмент, содержащий простую структуру, такую, как часть некоторого объекта и фон. Для того чтобы уменьшить влияние шума на наши наблюдения, следует выбрать ту область изображения, которая содержит полезный сигнал большой амплитуды. Используя яркости объекта и фона, мы приблизительно можем построить неразмытое изображение тех же размеров и с теми же особенностями, что и рассматриваемая часть исходного изображения. Обозначим рассматриваемую часть изображения как gs(x, y) и построенное изображение (которое в действительности представляет собой приближение для части неискаженного изображения

в рассматриваемой области) как f (x, y). Затем, имея в виду, что влияние шума пренебрежимо мало вследствие выбора области с большим полезным сигналом, на основе формулы (3) имеем [1]

H, (и, v) = .

Fs (и, v)

Далее, исходя из свойств функции Hs(u, v), можно сделать выводы о свойствах оптической передаточной функции системы H(u, v), принимая во внимание тот факт, что искажения в системе предполагаются транс-ляционно-инвариантными. Данный анализ поможет в оценке H(u, v).

а б

Рис. 4. Оценка искажающей функции с помощью импульса: а — яркий точечный световой импульс (показан с увеличением); б — реальное (искаженное) изображение светового импульса на выходе датчика ДЗЗ

Экспериментальный анализ. Можно провести точную оценку искажающей функции оптико-электронной системы, проведя наземные эксперименты с системой, аналогичной той, с помощью которой ведется съемка земной поверхности из космоса. Сначала необходимо подобрать параметры системы таким образом, чтобы искажения, получаемые в процессе эксперимента при формировании выходного сигнала в системе, как можно больше соответствовали искажениям, получаемым в данной системе при реальной космической съемке. Далее идея состоит в том, чтобы определить искажающую функцию путем получения изображения импульса (яркой точки очень малого размера) с использованием системы с подобранными значениями параметров. Импульс симулируется яркой световой точкой. Чтобы уменьшить влияние шума, яркость должна быть как можно больше и должны отсутствовать любые другие источники света (для этих целей можно использовать, например, лазерный луч, установленный в темной комнате). Направляя луч в объектив оптической системы, получаем некоторое размытое изображение. Затем, учитывая, что Фурье-преобразование импульса есть константа и пренебрегая воздействием шума, из формулы (3) получаем

О(и, V) Н(и, V) = 4 7, А

где G(u, V) — Фурье-преобразование полученного изображения; А — константа (Фурье-образ импульса), описывающая величину яркости светового импульса (соответствующий пример приведен на рис. 4).

Математическое моделирование. Математическое моделирование искажающей функции позволяет осознать саму суть процесса появления искажений, что необходимо для понимания работы алгоритмов восстановления изображений. Датчик дистанционного зондирования

Земли — сложная электронно-оптическая система, состоящая из целого ряда подсистем, таких как оптическая система формирования изображения, матрица светочувствительных сенсоров, аналого-цифровой преобразователь и др. Каждая из этих подсистем вносит свой вклад в общую искажающую функцию системы, что значительно усложняет математическое моделирование этой функции. Кроме того, имеют место и такие искажающие факторы, как бег подстилающей поверхности и пространственная ориентация платформы. Математическое моделирование искажающей функции помогает в борьбе с перечисленными факторами в заранее заданных пределах.

Рассмотрим математическое моделирование бега подстилающей поверхности, когда возникает смазывание изображения в результате равномерного поступательного движения изображения сцены относительно регистрирующей системы в процессе съемки. Предположим, что изображение f (х, у) движется, а функции х0(^), у0(/) определяют закон движения изображения по осям Ох и Оу соответственно. Полная экспозиция в любой точке записывающего носителя (например, пленки или матрицы сенсоров) определяется как интеграл по времени, т. е. времени, в течение которого открыт затвор регистрирующей системы, от величины мгновенной экспозиции.

Предположим, что затвор системы открывается и закрывается мгновенно и что процесс регистрации изображения (за исключением эффектов, связанных с движением) является идеальным. Тогда, если Т — время экспозиции (время интегрирования сигнала),

g(У) = {Дх - хоО), У - Уо 0)) Ж,

где g(x, у) — смазанное изображение.

Выполнив Фурье-преобразование функции от двух переменных g(x, у), получаем

0(и, V) = | |g(х, у)е-2п(ых+ууМхйу

Л

~т

{Дх - х0(У - Уо(0)с

_ о

Изменив порядок интегрирования, получим

т

-12 п( ыx+v

У) йхйУ.

О (ы, V) = j■

о

\ | /(х - хо(X), У - Уо(0) е~'2п(ых+"у>йхйУ

сИ. (4)

Выражение внутри скобок представляет собой Фурье-образ сдвинутой функции f (х - x0(t), y - y0(t)). Поэтому перепишем формулу (4) в следующем виде:

T T

G(u, v) = j>(u, v)e-2n(ux°(t}+vMt))dt = F(u, v)$e-2n(uxo(t}+vyo(t)]dt, (5) 0 0

причем второе равенство имеет место вследствие того, что функция F(u, v) не зависит от переменной t. Если пренебречь влиянием шума, то

G(u, v) = F(u, v)H(u, v). (6)

Тогда из формул (5) и (6) следует, что [1]

T

H(u, v) = je-2п(uxo(t}+vyo(t))dt.

0

Если функции x0(t), y0(t), определяющие закон движения изображения, известны, то передаточная функция H(u, v) может быть получена прямо из предыдущей формулы [1].

Алгоритмы восстановления изображений. Рассмотрим алгоритмы восстановления изображений, искаженных оператором H, пространственный образ которого (искажающая функция h) задан или определен с помощью методов, приведенных выше.

Инверсная фильтрация. Инверсная фильтрация является простейшим способом восстановления, который предполагает получение оценки F(u, v) Фурье-преобразования исходного изображения делением Фурье-образа искаженного изображения на оператор H:

F (u, v) = ^uA. (7)

H (u, v)

Деление в формуле (7) понимается как поэлементное. Подставив в формулу (7) выражение для G (u, v) из формулы (3), получим

F (u, v) = F (u, v) + Н(Ц. (8)

H(u, v)

Последнее выражение показывает, что, даже зная искажающую функцию, невозможно точно восстановить изображение. При использо-

вании данного метода неизбежно будет иметь место расхождение полученного результата с идеальным изображением, обусловленное слагае-

N(ц, V)

мым-за исключением случая, когда на искаженном изображении

Н (ы, V)

отсутствует аддитивный шум (Ы(и, V) = 0), что практически недостижимо при реальном восстановлении. Функция N (и, V) является Фурье-образом случайной величины (шума) и неизвестна. Кроме того, если функция Н (и, V) принимает нулевые или близкие к нулевым значения, N (ы, V)

то слагаемое-может на порядок превосходить значение функции

Н(ы, V)

¥(и, V) в данной точке частотной области, что приводит к значительным ошибочным результатам при восстановлении изображений данным методом. Иными словами, при таких условиях решение является неустойчивым [5]. Данный эффект часто встречается на практике, особенно при значительных показаниях шума. В этом состоит главный недостаток метода инверсной фильтрации. К его достоинствам можно отнести простоту реализации и низкую вычислительную сложность по сравнению с другими алгоритмами восстановления изображений.

Фильтрация методом минимизации среднеквадратического отклонения (фильтрация Винера). Данный метод восстановления изображений соединяет в себе учет свойств искажающей функции и статистических свойств шума. Фильтрация Винера основана на рассмотрении изображения и шума как случайных процессов. При этом задача

алгоритма состоит в том, чтобы найти такую оценку / для идеального (неискаженного) изображения / при которой среднеквадратическое отклонение этих величин было минимальным. Среднеквадратическое отклонение о задается формулой [1]:

а2 =

М {(/ - Л2},

где М{х} — функция математического ожидания аргумента х. Предполагается, что выполнены следующие условия [1]:

• шум и неискаженное изображение не коррелированы между собой;

• либо шум, либо неискаженное изображение имеют нулевое среднее значение;

• оценка линейно зависит от искаженного изображения.

При выполнении перечисленных условий минимум среднеквадра-тического отклонения достигается на функции, задаваемой в частотной области выражением [5]

Р (и, V) =

Н*(и, V) 8/(и, V)

8/ (и, V) Н(и, V)) + 8Л (и, V)

Н *(и, V) ^

О (и, V) =

Н(и, V) +

2 , (и,

Бг(и, V)

О (и, V) =

Н (и, V)!'

Н(и, V) , ,2 , 8Л (и, V)

Н (и,V)2 +

(и, V)

О (и, V),

(9)

где Бц(и, V) = |N(4, V) | — энергетический спектр шума; (и, V) = = |^(и, v)2| — энергетический спектр неискаженного изображения. В формуле (9) последнее равенство имеет место в силу того, что произведение комплексной функции на комплексно-сопряженную равно квадрату модуля этой функции. Данный результат был получен Н. Винером.

Анализируя выражение (9), можно заметить следующее:

• при отсутствии шума фильтр Винера переходит в обычный инверсный фильтр, рассмотренный ранее;

• при уменьшении спектральной плотности мощности исходного изображения передаточная функция фильтра Винера стремится к нулю;

• на частотах, соответствующих нулям передаточной функции формирующей системы, передаточная функция фильтра Винера также равна нулю. Таким образом, решается проблема сингулярности восстанавливающего фильтра. Отметим, что эта проблема присутствовала в алгоритме инверсной фильтрации (см. формулу (8)).

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

В случаях, когда спектры шума и неискаженного изображения неизвестны и не могут быть оценены, часто используется подход, состоящий в аппроксимации формулы (9) следующим выражением [1]:

Р(и, V) =

1

Н (и, V)!'

О (и, V),

Н(и, V) Н(и, V)!2 + К где К — определенная константа, соответствующая значению

(u, ^ 8/ (и, V)

в формуле (9). Она может быть подобрана в интерактивном режиме. Фильтрация, основанная на применении данного приема, называется параметрической фильтрацией Винера [1].

Фильтрация методом минимизации сглаживающего функционала со связью (метод Тихонова). Применение фильтрации Винера связано с дополнительной трудностью, состоящей в том, что энергетические спектры неискаженного изображения и шума также должны быть известны. Использование приближения (9) иногда позволяет получать хорошие результаты, однако применение константы в качестве оценки для отношения энергетических спектров не всегда приводит к удовлетворительному решению задачи.

Реализация метода минимизации сглаживающего функционала со связью требует только знания среднего значения и дисперсии шума. Это является важным преимуществом метода, поскольку обычно можно оценить указанные величины на основе заданного искаженного изображения [1].

Пользуясь определением свертки (2), можно записать выражение (1) в векторном виде:

Пусть изображение g(х, у) имеет размеры М х N. Тогда вектор g формируется таким образом, что первые М его элементов в точности совпадают с первой строкой изображения g(x, у), вторые М элементов — со второй строкой и т. д. Таким образом, в данном случае вектор g будет иметь длину МЫ. Такую же длину будут иметь векторы { и п. При этом матрица Н будет иметь размер МЫ х МЫ. Ее элементы определяются значениями h в свертке (1).

Формулировка задачи восстановления в матричном виде облегчает построение методов восстановления. Главной проблемой является чувствительность матрицы к шуму. Один из способов решения этой проблемы состоит в регуляризации задачи, которая заключается в замене исходной задачи на задачу нахождения экстремума некоторого сглаживающего функционала (данный подход также называется регуляризацией Тихонова). В качестве такого функционала можно использовать квадрат нормы лапласиана [1]

g = Н + п.

(10)

м-ш-1

с (/) (х, У ))2

2

(11)

х=о У=о

с дополнительным ограничением (связью) вида

22 § - ИТ =Н 2.

(12)

Решение оптимизационной задачи (11) в частотной области может быть представлено следующим выражением [6]:

( * Л р (и, V) = --Н^-¡2 О (и, V), (13)

^Н(и, V)!2 +у|Р(и, V)|2)

где параметр регуляризации у должен быть выбран в соответствии с (12), а функция P (и, V) есть Фурье-преобразование функции

Р (x, У ) =

0 -1 0" -1 4 -1 0 -10

т. е. функции, с помощью которой определяется оператор Лапласа [1].

Подчеркнем, что при обращении параметра регуляризации у в ноль выражение (13) сводится к инверсной фильтрации.

Значения параметра регуляризации у можно перебирать в интерактивном режиме до тех пор, пока не будет получен приемлемый результат. Однако для получения оптимального решения значение у должно быть выбрано таким образом, чтобы выполнялось условие связи (12). Ниже приведена итеративная процедура такого выбора [1].

Определим вектор невязки г следующим образом:

г = g - Ш. (14)

Поскольку решение Р(и, V) (и соответствующий ему вектор $), определяемое по формуле (13), есть функция от параметра у, вектор невязки г также зависит от этого параметра. Можно показать, что функционал невязки

ф( У) = гТ г = | И12

является монотонной функцией параметра у [1]. Нам необходимо выбрать параметр у таким образом, чтобы выполнялось условие

|| П II2 - а < II ^ II2 < II П II2 + а, (15)

где коэффициент а задает приемлемую точность выполнения результата связи. Если а = 0, имеет место II г Ц2 = II п II2 и вследствие (14) условие (12) выполняется точно.

Поскольку функционал невязки ф(у) является монотонной функцией параметра регуляризации, нахождение искомого значения у не представляет трудности. Один из алгоритмов состоит в следующем [7].

1. Задать начальное значение параметра регуляризации у.

2. Вычислить || г ||2.

3. Если условие (15) выполняется, цель достигнута, выходим из цикла. В противном случае необходимо увеличить значение у, если II г||2 < || П 1|2 - а, либо уменьшить его, если || г||2 > ||п ||2 + а. Повторить шаг 2, используя новую оценку у.

Для увеличения скорости сходимости можно использовать более продвинутые методы, например метод касательных Ньютона [1].

Для использования данного алгоритма необходимы величины II г ||2 и || п ||2. Для вычисления || г ||2 обратимся к формуле (14) и запишем ее в следующем виде:

Я (и, V) = G (и, V) - Н (и, V) Р( и, V),

причем функция г (х, у) может быть получена вычислением обратного Фурье-преобразования функции R(u, V). Итак,

М-1N-1

ИI2 = х, у ). (16)

х=о У=о

Рассмотрим дисперсию шума полного изображения

Л м-1 N-1

= ш Ц(П(^У)-тл)2, (П)

х=о У=о

где

М-1 N-1

тп = х, У) (18)

Ш1}/ х=о У=о

есть среднее значение шума.

Сопоставляя выражения (17), (18) с выражением для || п ||2, которое имеет тот же вид, что и (16), видим, что [8]

Н2 = МN (стП + т\).

Это важный результат, который показывает, что для реализации оптимального алгоритма восстановления методом минимизации сглаживающего функционала со связью достаточно знать лишь среднее значение шума и его дисперсию.

Эксперименты по восстановлению космических снимков. Проиллюстрируем работу рассмотренных алгоритмов восстановления изображений на примере реальных космических снимков. Для этого в качестве исходного изображения возьмем фрагмент снимка, полученного с космического аппарата Landsat-5 (рис. 5).

Рис. 5. Исходный снимок

Выполним программное размытие и зашумление исходного снимка. Результаты выполнения данных операций представлены на рис. 6 (для зашумления используется белый гауссов шум, дисперсия шума для всех

Рис. 6. Результаты программного размытия и зашумления исходного изображения:

а — размытого; б — размытого и слабозашумленного; в — размытого и среднеза-шумленного; г — размытого и сильнозашумленного

четырех снимков D = 0, математическое ожидание М для каждого снимка различно: для рис. 6, а M = 10-8, для 6, б M = 10-8, для 6, в M = 10-8 и для 6, г M = 10-8). Программа, выполняющая данное преобразование, реализована в среде MatLab. Размытие и зашумление осуществляется следующими командами:

PSF = fspecial('guassian', 15);

Blurred1=imnoise(Blurred, 'gaussian', noise mean, noise vari).

Метод инверсной фильтрации. Выполним восстановление размытых и зашумленных изображений методом инверсной фильтрации. Для восстановления изображения в разработанной программе используется процедура

imshow(deconvwnr(Blurred0, PSF, 0)); title( 'Результат восстановления размытого изображения').

Восстановленные методом инверсной фильтрации изображения представлены на рис. 7.

Рис. 7. Результаты восстановления изображения методом инверсной

фильтрации:

а — размытого; б — размытого и слабозашумленного; в — размытого и среднеза-шумленного; г — размытого и сильнозашумленного

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

Метод минимизации среднеквадратического отклонения (фильтрация Винера). Для оценки эффективности данного метода восстановления изображений выполним программное размытие и зашумление исходного изображения (см. рис. 5) со следующими параметрами белого гауссова шума: математическое ожидание М = 0, дисперсия Б = 10 5 (что на два порядка превышает уровень дисперсии шума, при котором инверсная фильтрация показала неприемлемые результаты восстановления — см. рис. 6, г). Результаты размытия и зашумления исходного изображения показаны на рис. 8 и 9.

Рис. 8. Размытое и зашумленное изображение

Рис. 9. Результат восстановления размытого и зашумленного изображения фильтрацией Винера (дисперсия шума

Проведем восстановление изображения методом фильтрации Винера. Для этого выполним следующую команду в разработанной программе:

deconvwnr(Blurred, PSF, estimated nsr).

Результат восстановления представлен на рис. 9. Очевидно, что данный алгоритм восстановления изображений обладает гораздо большей устойчивостью к шуму, чем инверсная фильтрация. Несмотря на более высокий показатель дисперсии шума (10-5) по сравнению с экспериментами с инверсной фильтрацией, алгоритм показал достаточно хорошие результаты восстановления изображения, в то время как восстановление методом инверсной фильтрации уже при уровне дисперсии шума 10-6 дало неприемлемые результаты.

Метод минимизации сглаживающего функционала со связью (регуляризация Тихонова). Проведем эксперимент по восстановлению размытого зашумленного изображения методом минимизации сглаживающего функционала со связью (регуляризация Тихонова). В качестве входных данных возьмем прежний снимок ДЗЗ и осуществим его размытие и зашумление программными средствами с точно такими же параметрами, как и в случае с фильтрацией Винера. После выполнения этих преобразований мы получим изображение, представленное на рис. 8. Затем выполним восстановление изображения методом регуляризации Тихонова. Программа, необходимая для выполнения этих действий, разработана в среде MatLab и приведена ниже.

I = im2double(imread('Fragment2.tif'));

figure(1); imshow(I); ^^е(ЛИсходное изображение');

PSF = fspecial('guassian', 15);

Blurred = imfilter(I, PSF, 'circular' , 'conv' );

noise mean = 0;

noise_var = 0.00001;

Blurred = imnoise(Blurred, 'gaussian', noise mean, noise var);

figure(2); imshow(Blurred); title('Pa3MbiToe и зашумленное изображение' );

estimated nsr = noise var / var(Blurred(:)); figure(4); imshow(deconvreg(Blurred, PSF, estimated nsr, [6e-4,1e10])); title('Tikhonov');

Результат выполнения этой программы (изображение, восстановленное методом регуляризации Тихонова) представлен на рис. 10.

Проведем визуальный анализ результатов восстановления расфокусированного и зашумленного изображения, представленного на рис. 9. Результаты восстановления, полученные при использовании

Рис. 10. Результат восстановления размытого и зашумленного изображения (дисперсия шума 10-5) с помощью регуляризации Тихонова

фильтрации методом минимизации среднеквадратического отклонения (фильтрация Винера) и фильтрации методом минимизации сглаживающего функционала со связью (регуляризация Тихонова), представлены на рис. 8 и 9 соответственно. Можно отметить, что результаты восстановления изображений этими двумя методами сопоставимы один с другим, но в данном случае регуляризация Тихонова дала более четкий и приближенный к оригиналу результат.

Выводы

1. Для осуществления радиометрической коррекции космических снимков необходимо обладать информацией об искажающей функции датчика ДЗЗ, для построения которой требуется знание параметров конкретного датчика (разрешение и размер светочувствительной матрицы, фокусное расстояние линзы, временной интервал дискретизации и др.), а также параметров орбиты космического аппарата. По причине отсутствия информации о параметрах конкретных датчиков ДЗЗ в работе использовалось программное зашумление, размытие и смазывание изображения с его последующим восстановлением различными методами.

2. Алгоритм восстановления изображений методом инверсной фильтрации показал идеальные результаты при восстановлении размытых незашумленных изображений, однако даже при небольшом уровне шума (гауссов шум с математическим ожиданием 0 и дисперсией 10 8) на восстановленном изображении стали проявляться значительные искажения. При уровне дисперсии шума 10 6 результат восстановления стал практически неразборчивым, что говорит о низкой стойкости этого метода к аддитивному шуму.

3. При значении дисперсии аддитивного шума 10 5 алгоритм восстановления изображений методом параметрической фильтрации Ви-

нера показал относительно хорошие результаты. Это свидетельствует о гораздо более высокой стойкости данного метода к аддитивному шуму по сравнению с инверсной фильтрацией.

4. Наиболее приближенный к оригиналу результат восстановления размытого зашумленного изображения (дисперсия шума 10-5) показал алгоритм восстановления изображений методом минимизации сглаживающего функционала (регуляризация Тихонова).

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

Работа выполнена при финансовой поддержке Министерства образования и науки Российской Федерации в рамках мероприятия 1.9 Федеральной целевой программы «Исследования и разработки по приоритетным направлениям развития научно-технического комплекса России на 2007-2013 годы», государственный контракт от 18 августа 2011 г. № 11.519.11.5002, шифр 2011-1.9-519-004.

ЛИТЕРАТУРА

[1] Гонсалес Р., Вудс Р. Цифровая обработка изображений, Москва, Техносфера, 2005, 1072 с.

[2] Разработка методического и алгоритмического обеспечения проведения мониторинга качественных и количественных характеристик локальных экосистем с использованием современных и перспективных космических средств дистанционного зондирования Земли. Отчет о НИР Минобрнауки РФ «Исследование возможностей современных и перспективных космических средств дистанционного зондирования Земли для мониторинга состояния локальных экосистем» в рамках мероприятия 1.9 Федеральной целевой программы «Исследования и разработки по приоритетным направлениям развития научно-технического комплекса России на 2007-2013 годы».

[3] Шовенгердт Р.А. Дистанционное зондирование. Модели и методы обработки изображений. Москва, Техносфера, 2010, 560 с.

[4] Численное моделирование мониторинга качественных и количественных характеристик локальных экосистем с использованием современных и перспективных средств дистанционного зондирования Земли по разработанной методике и анализу результатов численного моделирования. Отчет о НИР Минобрнауки РФ «Исследование возможностей современных и перспективных космических средств дистанционного зондирования Земли для мониторинга состояния локальных экосистем» в рамках мероприятия 1.9 Федеральной целевой программы «Исследования и разработки по приоритетным направлениям развития научно-технического комплекса России на 2007-2013 годы».

[5] Ягола А.Г. Некорректные задачи и методы их численного решения. Спец. курс для аспирантов МГУ им. М.В. Ломоносова, Москва, 2005, 21 с.

[6] Гонсалес Р., Вудс Р., Эддинс С. Цифровая обработка изображений в среде MathLab. Москва, Техносфера, 2006, 618 с.

[7] Винер Н. Нелинейные задачи в теории случайных процессов. Москва, Иностранная литература, 1961, 158 с.

[8] Петров Ю.П., Сизиков В.С. Корректные, некорректные и промежуточные задачи с приложениями. Санкт-Петербург, Политехника, 2003, 261 с.

Статья поступила в редакцию 21.05.2013

Ссылку на эту статью просим оформлять следующим образом: В.И. Майорова, А.М. Банников, К.И. Зайцев. Математическое моделирование процесса радиометрической коррекции снимков ДЗЗ. Инженерный журнал: наука и инновации, 2013, вып. 3. URL: http://engjoumal.ru/catalog/mathmodel/hidden/641.html

Майорова Вера Ивановна родилась в 1956 г., окончила МГТУ им. Н.Э. Баумана в 1979 г. Д-р техн. наук, проф. кафедры «Космические аппараты и ракеты-носители» МГТУ им. Н.Э. Баумана. Автор более 70 работ в области разработки космической техники и образовательных космических технологий. e-mail: victoria.mayorova@ gmail.com

Банников Алексей Михайлович родился в 1991 г Студент 5-го курса кафедры «Информационная безопасность» МГТУ им. Н.Э. Баумана. Автор двух публикаций в области применения информационных технологий в космонавтике. e-mail: alexm.fighter@ gmail.com

Зайцев Кирилл Игоревич родился в 1989 г., окончил МГТУ им. Н.Э. Баумана в 2012 г. Аспирант кафедры «Лазерные и оптико-электронные системы» (РЛ-2) факультета «Радиоэлектроника и лазерная техника» МГТУ им. Н.Э. Баумана. Автор 15 публикаций. Сфера научных интересов: цифровая обработка изображений; стереоскопические методы восстановления трехмерных образов объектов; использование терагерцовой спектроскопии и терагерцовых изображающих систем в целях медицинской диагностики; обратные задачи оптики; вычислительная электродинамика. e-mail: kirzay@gmail.com