Обработка и фильтрация данных дистанционного зондирования

1.3 Атмосферная коррекция

1.4 Восстановление пропущенных пикселов

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

Спутники функционируют на орбите в течение нескольких лет, с течением времени измерительная аппаратура деградирует под воздействием неблагоприятных факторов космического пространства. Поэтому показания датчиков сканеров необходимо корректировать. Эта процедура носит название радиометрической коррекции. Рассмотрим радиометрическую коррекцию на примере обработки данных сканера AVHRR спутника NOAA1.

Сканер AVHRR имеет по одному каналу видимого и ближнего инфракрасного диапазона (1-й и 2-й каналы), один ИК-канал на 3,5 мкм (3-й канал) и два канала теплового ИК-диапазона (4-й и 5-й каналы). Предусмотрена бортовая калибровка последних трех каналов, для этого сканер направляется на космическое пространство (эта точка принимается за ноль) и в полость абсолютно черного тела, установленного на борту, что дает две точки для коррекции температуры по линейному закону. Корректировочные коэффициенты для 3-5-го каналов включены в файл данных, передаваемых со спутника; 1-й и 2-й каналы калибруются только на Земле перед запуском спутника, соответствующие коэффициенты также включены в файл данных. С течением времени возникает необходимость корректировать показания 1-го и 2-го каналов, применяется процедура, состоящая в том, что наблюдаются одни и те же объекты на Земле, текущая интенсивность излучения сравнивается с результатами предыдущих наблюдений. И текущие, и полученные прежде данные подвержены разбросу из-за влияния атмосферы, прозрачность которой постоянно меняется, из-за естественной изменчивости природных объектов и других факторов. Корректировочные коэффициенты для 1-го и 2-го каналов определяются путем статистической обработки. Коэффициенты ежемесячно обновляются, их можно найти в сети Интернет2.

В качестве примера рассмотрим процедуру калибровки и корректировки данных 4-го канала. Показания U датчика этого канала после аналого-цифрового преобразования лежат в пределах 0-1 023 (10-битное кодирование). Пусть U = 427. Плотность потока мощности можно найти по формуле

где A = 590,888 Вт-мкм/(м2*ср),

K = -1,60156 Вт-мкм/(м2*ср) - корректировочные коэффициенты из файла данных.

В нашем примере B(л, T) = 907,022/л2 Вт/(м2*мкм*ср).

Из формулы Планка3 получаем, что радиационная температура при л=10,96 мкм, с1=1,1911*108 Вт*мкм /(м2*ср), c2=14 388 мкм*К

T = л/c2 ln (c1/ л5B + 1) = 284,57 К.

Пусть U = 428, т.е. показания датчика изменились на единицу. В этом

случае В = 905,421/л2 Вт/(м*мкм*ср), Т = 284,46 К Таким образом, изменению показаний датчика ДU = 1 соответствует изменение температуры ДT = 0,1 К. Величина ДT характеризует радиационное разрешение датчика.

1.3 Атмосферная коррекция

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

Лучшим выходом из положения была бы установка по всей поверхности суши обширной сети солнечных спектрофотометров. Назначение этих приборов - измерение в различных участках спектра интенсивности солнечного излучения I, прошедшего через атмосферу. Зная интенсивность Iо за пределами атмосферы, по закону Бугера

I = Io exp(-ф sec д)

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

ф = 1/sec д * lnI0/I.

Здесь д - зенитный угол Солнца.

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

Иначе дело обстоит с атмосферой над морями и океанами. В красном и ИК-участке спектра поверхность воды по своим оптическим характеристикам близка к абсолютно черному телу. Существенно больший коэффициент отражения и рассеяния имеют туманы, дымки, облака, их хорошо видно на фоне воды. Это позволяет оценить оптическую толщину ф. Данные по ф над морями и океанами помещены в сети Интернет4, их можно использовать для коррекции спутниковых изображений прибрежных районов.

На рис. 4, а приведено изображение части Хабаровского края (слева), Татарского пролива и о. Сахалин (справа), полученное 14 августа 1998 г. (1-й канал сканера AVHRR спутника NOAA).

А Б

Рис. 4. Атмосферная коррекция изображения прибрежных районов

Это изображение очень низкого качества. В это время на территории Хабаровского края и на Сахалине бушевали обширные лесные пожары, весь район сильно задымлен. Была произведена атмосферная коррекция с использованием закона Бугера, данные об оптической толщине ф взяты из сети Интернет, ф лежало в пределах 0,03-0,23, угол д изменялся от минус 8 до 14°. Результаты коррекции показаны на рис. 4, б. После коррекции просматривается Татарский пролив и несколько лучше (чем на рис. 4, а) - о. Сахалин. Для Хабаровского края использованы данные ф по Татарскому проливу и Охотскому морю, результат, как и следовало ожидать, оказался неудовлетворительным: в левом верхнем углу появилась темная область, похожая на поверхность моря, хотя на самом деле это суша.

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

Алгоритм порогового выделения загрязненных областей является наиболее простым, но достаточно эффективным. Идея алгоритма основана на выделении границы между «загрязненной» и «чистой» областью, после чего все значения изображения в пределах этой границы исключаются из рассмотрения. Порог для 1-го и 2-го каналов сканера AVHRR спутника NOAA можно выбрать, если вычислить значения нормализованного вегетационного индекса NDVI для всего изображения. Эта процедура основана на чувствительности NDVI к присутствию облаков и аэрозолей на изображении. Для растений в нормальном состоянии NDVI близок к 0,6-0,65, для растений в угнетенном состоянии - к 0,3 - 0,4. Облака и аэрозоли в 1-м и 2-м каналах имеют приблизительно одинаковую яркость, поэтому NDVI близок к 0. Таким образом, пороговыми являются низкие значения NDVI. Порог для 3, 4 и 5-го каналов выбирается как среднее значение температуры поверхности и облаков.

Кроме солнечного излучения, отраженного и рассеянного поверхностью, сканеры спутника принимают солнечное излучение, рассеянное на молекулах газов и рассеянное на аэрозолях. Интенсивность молекулярного рассеяния обратно пропорциональна длине волны в четвертой степени, эффект наиболее заметен в коротковолновой части спектра, он ответствен за голубой цвет неба. Рассеяние на аэрозолях (размер частиц от 0,1 л до 10л) приводит к зависимости интенсивности от длины волны л-б, 0 < б < 4. Частицы дыма, пыли, облаков имеют размеры много большие, чем длина волны видимого и ИК-диапазонов, здесь нет зависимости интенсивности рассеянного солнечного излучения от длины волны. Рассеяние в атмосфере приводит к дополнительной засветке, заметной в 1-м и 2-м каналах сканера AVHRR спутника NOAA и еще более заметной в каналах сканеров спутников «Ресурс-О», LANDSAT, SPOT, работающих в зеленой области видимого диапазона. Засветка различна для разных каналов, она способна привести, например, к ошибкам в оценке вегетационного индекса растительного покрова.

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

Рис. 5. Гистограмма уровней яркости при наличии рассеяния в атмосфере

Методика коррекции состоит в том, что из всех уровней яркости вычитается величина ДI, определяемая по гистограмме.

1.4 Восстановление пропущенных пикселов

Пусть, например, уровни некоторого черно-белого изображения занимают интервал от 6 до 158 со средним значением яркости 67 при возможном наибольшем интервале значений от 0 до 255. На рис. 7, а приведена гистограмма яркостей исходного изображения, показывающая, сколько пикселов N с близким значением яркости f попадает в интервал от fi до f +Дf

Рис. 7. Исходное изображение (а) и изображение после линейной растяжки гистограммы (б)

Рис. 8. Нормализация гистограммы

Это изображение является малоконтрастным, превалирует темный оттенок. Возможным методом улучшения контраста может стать так называемая линейная растяжка гистограммы (stretch), когда уровням исходного изображения, лежащим в интервале [fmin, fmax], присваиваются новые значения с тем, чтобы охватить весь возможный интервал изменения яркости, в данном случае [0, 255]. При этом контраст существенно увеличивается (рис. 7, б). Преобразование уровней яркости осуществляется по формуле

gi = с + dfi,

Выражение 4.

где fi - старое значение яркости i-го пиксела; gi - новое значение; c, d - коэффициенты. Для рис. 7, а fmin = 6,fmax = 158. Выберем c и d таким образом, чтобы gmin = 0, gmax = 255. Из формулы 4 получаем c = -10,01; d = 1,67.

Еще более можно улучшить контраст, используя нормализацию гистограммы (рис. 8), когда на весь максимальный интервал уровней яркости [0, 255] растягивается не вся гистограмма, лежащая в пределах от fmin до fmax, а её наиболее интенсивный участок (fmin', fmax'), из рассмотрения исключаются малоинформативные «хвосты». На рис. 8, б исключено 5 % пикселов.

Целью выравнивания гистограммы (эту процедуру называют также линеаризацией и эквализацией - equalization) является такое преобразование, при котором в идеале все уровни яркости приобрели бы одинаковую частоту, а гистограмма яркостей отвечала бы равномерному закону распределения (рис.9).

Рис. 9. Гистограмма, отвечающая равномерному закону распределения

Пусть изображение имеет формат: N пикселов по горизонтали и M по вертикали, число уровней квантования яркости равно J. Общее число пикселов равно N * M, на один уровень яркости попадает в среднем n0 = N * M/J пикселов. Например, N = M = 512, J = 256. В этом случае n0 = 1 024. Расстояние Дf между дискретными уровнями яркости от fi до f i+1 в гистограмме исходного изображения одинаковое, но на каждый уровень выпадает различное число пикселов. При эквализации гистограммы расстояние Дgi между уровнями gi и gi+1 различно, но число пикселов на каждом уровне в среднем одинаковое и равно no. Алгоритм эквализации несложен. Пусть уровнями с малой яркостью обладает небольшое количество пикселов, как на рис. 10, а. Например, уровень яркости 0 на исходном изображении имеют 188 пикселов, уровень 1 - 347 пикселов, уровень 2 - 544 пиксела. В сумме это 1 079 пикселов, т. е. приблизительно n0. Присвоим всем этим пикселам уровень 0. Пусть на исходном изображении число пикселов с уровнями яркости 3 и 4 в сумме приблизительно также равно n0. Этим пикселам присваивается уровень 1. С другой стороны, пусть число пикселов с уровнем 45 на исходном изображении составляет 3 012, т. е. приблизительно 3n0. Всем этим пикселам присваивается некоторый одинаковый уровень gi, не обязательно равный 45, а соседние два уровня остаются незаполненными. Эти процедуры повторяются для всех уровней яркости. Результат эквализации можно видеть на рис. 10, б.

а б

Рис. 10. Эквализация гистограммы

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

Процедуру видоизменения гистограммы можно рассматривать как по пиксельное преобразование входной яркости fj, f0 <= f <= fj, в выходную яркость gk, g0 <= gk <= gK, в результате которого исходное распределение вероятностей P{fj} переходит в распределение вероятностей P{gk}, имеющее желаемую форму. Очевидно, что сумма вероятностей яркостей всех пикселов должна равняться единице:

Вероятность попадания исходной яркости f i в интервал от 0 до m должна равняться вероятности попадания яркости преобразованного изображения gk в интервал от 0 до n для всех m <= J, n < =K:

В случае конкретного изображения распределение в левой части заменяют на гистограмму H(fj), поэтому

Решая это уравнение, можно найти требуемое преобразование gk=T{fj}. Решение записывается в виде таблицы, в которой для каждого входного уровня fj указывается соответствующий выходной уровень gk.

2. Линейная пространственно-инвариантная фильтрация изображений

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

Выражение 5

при h1(x, y) ~ exp [-(x2 + y2)/б].

2.1 Пространственные инвариантные операторы

Модель искаженного помехами непрерывного изображения имеет вид

f(x, y) = m(x, y) * Fs(x, y) + n(x, y),

где f(x, y) - искаженное изображение; m(x, y) - мультипликативная помеха, модулирующая изображение по яркости; F - функционал, описывающий геометрические и радиометрические искажения, а также расфокусировку; s(x, y) - исходное изображение; n(x, y) - аддитивная помеха, накладывающаяся на изображение. Модуляция спутникового изображения по яркости может наблюдаться из-за того, что атмосфера над различными точками Земли имеет разную прозрачность, восходящее излучение от этих точек проходит различный путь в атмосфере.

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

m(x, y) = 1. Таким образом,

f(x, y) = Fs(x, y) + n(x, y).

Результат реставрации s(x, y) = g(x, y) запишем как следствие воздействия на f(x, y) некоторого оператора:

g(x, y) = Tf(x, y).

Оператор T (системный оператор) указывает на правило, по которому «входному сигналу» f(x, y) ставится в соответствие «выходной сигнал» g(x, y). Для того чтобы модель была полной, необходимо также указать области допустимых значений f(x, y) и g(x, y). При реставрации применяют оператор Т, минимизирующей расстояние между g(x, y) и s(x, y) при заданных статистических характеристиках случайных полей s(x, y), n(x, y) и известном F. В качестве критерия близости g(x, y) и s(x, y) часто используют критерий минимума среднеквадратической ошибки:

min <[g(x, y) - s(x, y)]2 >.

В задачах улучшения изображений обычно считается, что n(x, y) = 0, функцией оператора Т является сглаживание резких перепадов яркости, подчеркивание или выделения контуров и т. п.

Мы будем рассматривать пространственно-инвариантные операторы, выходная реакция которых не зависит от изменения начала отсчета по x и по y и от ориентации объектов на изображении. Первое условие означает, что оператор переводит однородное случайное поле в однородное. Второе условие означает, что оператор переводит изотропное поле в изотропное. Отметим, что свойства пространственной инвариантности выполняются строго, если области допустимых значений координат x, y попадают в интервал (-?, ?). Реальные изображения имеют конечные размеры, A <= x <= B; C <= y <= D, условие пространственной инвариантности выполняется приближенно.

Оператор называется линейным, если для него справедлив принцип суперпозиции - реакция на сумму сигналов f1(x, y) и f2(x, y) равна сумме реакций на каждое из воздействий в отдельности

T (f1(x, y) + f2(x, y)) = T f1(x, y) + Т f2(x, y),

для любого произвольного числа б справедливо:

Tбf(x, y) = бTf(x, y).

Свойства линейности выполняются строго, если области допустимых значений яркости f, g попадают в интервал (-?, ?). При цифровой обработке яркость - величина вещественная, неотрицательная и ограниченная, обычно 0<=f, g<= 255. Если каждому g(x, y) отвечает единственное f(x, y), то оператор Т может быть представлен в виде функционального степенного ряда (ряда Вольтера).

Выражение 5

Здесь интегрирование ведется по всей области, где определены x, y; записаны два члена ряда Вольтера (линейный и квадратичный); весовые множители

называются ядрами Вольтера первого и второго порядка.

Рис. 11. Пример функции рассеяния точки

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

Ядро первого порядка h1(x, y, x', yr) в оптике именуют функцией рассеяния точки (ФРТ). Это изображение точечного источника на выходе оптической системы, которое уже является не точкой, а некоторым пятном. В соответствии с выражением 5, первый член, все точки изображения f(x', y') превращаются в пятна, происходит суммирование (интегрирование) всех пятен. Не следует думать, что эта процедура обязательно приводит к расфокусировке изображения, наоборот, можно подобрать такую ФРТ, которая позволит сфокусировать расфокусированное изображение. На рис. 11 представлена одна из возможных ФРТ.

Для того чтобы для ФРТ выполнялось условие пространственной инвариантности, т. е. чтобы ФРТ не изменялась при изменении начала отсчета по x и по y, она должна иметь вид h1(x, y, x', yr) = h1(x -x', y -y') В этом случае h1(x, y, x', y') = h1(x + x0, y + y0, x' + x0, y' + y0 ). Кроме того, ФРТ должна обладать осевой симметрией.

При обработке растровых изображений на прямоугольной сетке проще всего реализовать ФРТ конечных размеров в виде прямоугольной матрицы форматом NxN, например, 3x3:

Только три элемента матрицы размером 3x3 независимы, в этом случае матрица инвариантна относительно поворотов, кратных 90°. Опыт обработки изображений показывает, что отсутствие более строгой осевой симметрии ФРТ слабо сказывается на результатах. Иногда используют 8-угольные матрицы, инвариантные относительно поворотов на 45°.

2.2 Линейные преобразования в частотной плоскости

Реальное растровое изображение fn,m имеет конечные размеры A<= x<=B, C<=y<=D и состоит из отдельных пикселов, расположенных с некоторым шагом в узлах прямоугольной сетки. В этом случае для перехода в частотную плоскость применяется двойное дискретное преобразование Фурье (ДДПФ):

Здесь N - число пикселов по x; M - число пикселов по y. Величины Fp,qназываются коэффициентами ДДПФ, которые вычисляют следующим образом: вначале в выражении (12) производится суммирование по n (по строкам), потом полученный числовой массив суммируется по m (по столбцам), можно и наоборот. Разработан также алгоритм, сводящий двумерное дискретное преобразование Фурье к одномерному.

Коэффициент Fo,o - это среднее значение яркости изображения. Общее число комплексных коэффициентов Fp, q равно N * M, однако часть из них связана между собой. Если N - четное число, например, то

звездочка означает комплексную сопряженность, так что число независимых коэффициентов равно N * M/4.

Для вычисления ДДПФ согласно выражению (12) необходимо выполнить N2 * M2 операций с комплексными числами, что требует значительных затрат времени. Если N = M = 512, то общее число операций составляет 5124 ? 6,81010. С целью ускорения вычислений разработаны алгоритмы быстрого преобразования Фурье (БПФ), отличающиеся логарифмической сложностью. При их применении общее число операций оценивается как N logN * M logM, что для приведенного примера составляет 2*107 операций. Таким образом, вычисления ускоряются более чем в 3 000 раз.

Глобальная линейная фильтрация изображений в частотной плоскости осуществляется умножением Fp, q на частотный коэффициент передачи Кр, q. На рис. 12, а приведено изображение г. Красноярска (спутник LNDSAT, разрешение 30 м), на рис. 12, б - модуль ДДПФ этого изображения. На рис. 12, в показан результат глобальной фильтрации фильтром с коэффициентом передачи Кр, q = exp [-(р2 + q2)/r], который изображен на рис. 12, г (фильтром нижних частот). Результат обработки фильтром верхних частот с коэффициентом передачи Кр, q = 1 - exp [- (р + q )/r], который изображен на рис. 12, е, приведен на рис. 12, д.

Темный тон на изображении фильтров означает, что фильтр не пропускает помеченные этим тоном частоты, белый означает пропускание. Фильтр нижних частот сглаживает изображение, фильтр верхних частот выделяет контуры объектов на изображении.

Рис. 12. Глобальная фильтрация: а - исходное изображение; б - модуль ДДПФ этого изображения; в - результат обработки фильтром нижних частот; д - результат обработки фильтром верхних частот; г - фильтр нижних частот; е - фильтр верхних частот

2.3 Применение линейной фильтрации в частотной плоскости

При линейной фильтрации изображений в частотной плоскости требуется умножить спектр пространственных частот на частотный коэффициент передачи согласно выражению (11) и выполнить два двумерных преобразования Фурье с использованием алгоритма БПФ - прямое и обратное. Преобразование Фурье осуществляется от всего изображения целиком, спектр F(u, v) хранит информацию обо всем изображении целиком, поэтому мы говорим о глобальной фильтрации. В зависимости от выбора коэффициента передачи К(и,v) можно выделять изображение на фоне помех, улучшать и ухудшать резкость изображения, выделять контуры объектов на изображении и т. п. Качество обработки при этом несколько лучше, чем при локальной линейной фильтрации. Типичной задачей глобальной фильтрации, не характерной для цифровой обработки космических изображений земной поверхности, является восстановление расфокусированных изображений. Запишем спектр такого изображения в виде

F1(u, v) = F(u, v) * К(и, v),

где F(u, v) - спектр исходного изображения; К(и, v) - коэффициент передачи оптической системы, отвечающий ФРТ вида рис. 11. Примеры расфокусированных изображений приведены на рис. 13.

Рис. 13. Обработки перепада яркости маской Н2

Очевидно, что для того чтобы определить F(u, v) по известным F1(u, v) и К(и, v), необходимо умножить F1(u, v) на К1(и, v) = 1/К(и, v). Эта процедура носит название инверсной фильтрации. Задача инверсной фильтрации относится к числу обратных некорректных задач математической физики. Во-первых, решение может не существовать. Во-вторых, если решение существует, то может быть не единственным. Заметим, что функция К(и, v) в некоторых точках может быть равна нулю, тогда решение обращается в бесконечность. В-третьих, решение может быть неустойчивым, т. е. небольшие вариации исходных данных могут привести к существенным изменениям решения. Выходом из положения является регуляризация решения. Коэффициент передачи записывается в форме

К1(и, v) = К(и, v) / [К2(и, v) + б г(u, v)].

Выражение 14.

Здесь г(и, v) - стабилизирующая функция, например, г(и, v) = и2 + v2 ; б - параметр регуляризации. Подбором б обычно удается достаточно качественно восстановить расфокусированное изображение.

Фильтрация в частотной плоскости - эффективное средство для восстановления изображений, искаженных аддитивным шумом. Пусть изображение s(x, y) наблюдается на фоне некоррелированного с s(x, y) шума n(x, y):

f(x, y) = s(x, y) + n(x, y).

Предполагается, что f(x, y), s(x, y) и n(x, y) - однородные и изотропные гауссовские случайные поля с функциями корреляции Rf (r), Rs(r), Rn(r), соответственно, с нулевыми, средними и ограниченными дисперсиями, Rs(r)=<s(x, y)s(x1, y1)>, остальные функции корреляции определены аналогичным образом; аргумент r = v(x - x1 )2 +(y - y1 )2 . Из математической статистики известно, что если и полезный сигнал (в нашем случае это изображение), и шум имеют нормальный (гауссовский) закон распределения, то самый лучший эффект фильтрации, по крайней мере, с учетом минимума среднего квадрата ошибки обеспечивает линейный фильтр. Для выделения изображения на фоне шума используем линейное преобразование

Требуется установить ФРТ h1(x, y), обеспечивающую минимум среднего квадрата ошибки фильтрации е2 = <|g(x, y) - s(x, y)|2 >. Находим, что

где

учтено, что <s(x, y)n(x1, y1)> = 0.

Пусть h1 мало отличается от оптимального значения h1*, т. е. h1 = h1* + б-в(x, y), здесь б - неопределенный множитель Лагранжа; в(x, y) - некоторая функция. Необходимое условие минимума состоит в том, что dе/dб = 0 при б=0.

Таким образом,

Это равенство должно выполняться при произвольной в(x - x', y - y'), откуда

Выражение 15. Перейдем от функций корреляции к спектру мощности G(u, v) и от ФРТ к частотному коэффициенту передачи К(и, v). Спектр мощности связан с функцией корреляции преобразованием Фурье:

Выражение 16

Подстановка выражений (10) и (13) при Gf(u, v) = Gs(u, v) + Gn(u, v) в формулу (15) дает

оптимальный коэффициент передачи

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

В случае, когда изображение вначале искажено из-за воздействия некоторого линейного оператора, а затем на него наложен шум, возможно объединение функций фильтра Винера и инверсного фильтра. Коэффициент передачи равен

K2(u, v) = 1/[i + К12(u, v)/K*(u, v)],

где К1(u, v) определяется выражением (14) при б = 0. Гомоморфный фильтр имеет коэффициент передачи

K3(u, v) = K1(u, v)/ v1+ K21(u, v)/K *(u, v).

Его также называют обобщенным линейным фильтром.

2.4 Линейная локальная фильтрация

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

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

суммирование ведется по некоторой окрестности D точки (i, j); бkl -- значения ФРТ в этой окрестности. Яркости пикселов f в этой точке и в её окрестности умножаются на коэффициенты бkl, преобразованная яркость (i,j)-гo пиксела есть сумма этих произведений. Обычно набор коэффициентов бkl представляют в виде прямоугольной матрицы (маски), например, размерности 3x3:

Элементы матрицы удовлетворяют условию пространственной инвариантности, поэтому a11 = a13 = a31 = a33, a12 = a21 = a23 = a32 .

Фильтрация осуществляется перемещением слева направо (или сверху вниз) маски на один пиксел. При каждом положении апертуры производятся упомянутые выше операции, а именно перемножение весовых множителей бkl с соответствующими значениями яркостей исходного изображения и суммированием произведений. Полученное значение присваивается центральному (i,j)-му пикселу. Обычно это значение делится на заранее заданное число K (нормирующий множитель). Маска содержит нечетное число строк и столбцов, чтобы центральный элемент определялся однозначно.

Рассмотрим некоторые фильтры, сглаживающие шум. Пусть маска размером 3x3 имеет вид

Тогда яркость (i, j)-ro пиксела после фильтрации определяем по формуле

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

У фильтров Н1--Н4 нормирующие множители K подобраны таким образом, чтобы не происходило изменения средней яркости обработанного изображения. Наряду с масками 3x3 используются маски большей размерности, например, 5x5, 7x7 и т. п. В отличие от фильтра Н2, у фильтров Н1, Н3, Н4 весовые коэффициенты на пересечении главных диагоналей матрицы больше, чем коэффициенты, стоящие на периферии. Фильтры Н1, Н3, Н4 дают более плавное изменение яркости по изображению, чем Н2.

Пусть отсчеты полезного изображения мало меняются в пределах маски, а отсчеты аддитивного шума случайны и независимы или слабо зависимы со статистической точки зрения. В этом случае механизм подавления шума с использованием приведенных фильтров состоит в том, что при суммировании шумы компенсируют друг друга. Эта компенсация будет происходить тем успешнее, чем большее число членов в сумме, т. е. чем больше размер (апертура) маски. Пусть, например, используется маска NxN, в пределах её полезное изображение имеет постоянную яркость f, шум -- аддитивный, с независимыми значениями отсчетов nk,m, средним значением м = 0 и дисперсией у2 в пределах маски (такой шум называют белым). Отношение квадрата яркости (i, j)-гo пиксела к дисперсии шума, т.е. отношение сигнал/шум, равно f2 /у2 .

Рассмотрим, например, маску типа Н2:

Средний квадрат яркости равен f2, средний квадрат интенсивности шума

Двойная сумма отвечает k = p, m = q, эта сумма равна у2/N2 . Четырехкратная сумма равна нулю, так как отсчеты шума при k ? p, m ? q независимы: <nk,m np,q> = 0. Таким образом, в результате фильтрации отношение сигнал/шум становится равным N2f 2/у2 , т. е. возрастает пропорционально площади маски. Отношение яркости (i, j)-гo пиксела полезного изображения к среднеквадратическому отклонению шума возрастает пропорционально N2. Следовательно, использование маски 3x3 в среднем повышает отношение сигнала к шуму в 9 раз.

При импульсной помехе механизм подавления состоит в том, что импульс «расплывается» и становится мало заметным на общем фоне.

Однако часто в пределах апертуры значения полезного изображения все же изменяются заметным образом. Это бывает, в частности, когда в пределы маски попадают контуры. С физической точки зрения, все H1-H4 являются фильтрами нижних частот (усредняющими фильтрами), подавляющими высокочастотные гармоники и шума, и полезного изображения. Это приводит не только к ослаблению шума, но и к размыванию контуров на изображении. Пусть, например, на изображении имеется перепад яркости от 100 к 200 (рис. 14), изображение обрабатывается маской Н2.

Если маска находится в положении а, когда в её пределы попадают только значения яркости f = 100, то gi,j = 1/9 * 9 * 100 = 100. Если маска находится в положении г, то gi, j = 1/9 * 9 * 200 = 200. Таким образом, в этих случаях изменения яркости не происходит.