Содержание к диссертации
Введение
ГЛАВА 1. Анализ состояний исследований в области восстановления изображений земной поверхности 8
1.1. Анализ формирования спутниковых изображений земной поверхности 8
1.2. Проблемы обработки изображений земной поверхности при наличии функции рассеяния точки 17
1.3. Анализ существующих методов восстановления изображений земной поверхности, полученных в оптическом диапазоне 26
1.4. Постановка задач исследования 29
Выводы к главе 1 32
ГЛАВА 2. Разработка метода восстановления спутниковых изображений земной поверхности на основе модификации оператора тихонова 33
2.1. Аддитивное представление восстанавливающего оператора 33
2.2. Адаптивное оценивание параметров восстанавливающего оператора 42
2.3. Итерационная процедура восстановления изображений земной поверхности 49
2.4. Определение параметров восстанавливающего оператора для эталонных снимков участков земной поверхности 53
Выводы к главе 2 58
ГЛАВА 3. Разработка алгоритмов восстановления изображений земной поверхности 59
3.1. Алгоритм численной реализации метода восстановления 59
3.2. Исследование алгоритма восстановления на модельных примерах 65
3.3. Алгоритм восстановления изображений при наличии паразитных полос 73
3.4. Алгоритм восстановления при геометрическом несовершенстве фотоприёмных линеек 82
Выводы к главе 3 91
ГЛАВА 4. Программная реализация разработанных алгоритмов восстановления изображений земной поверхности 92
4.1. Разработка прототипа программной реализации алгоритмов восстановления изображений земной поверхности 92
4.2. Проведение сравнительных вычислительных экспериментов 103
Выводы к главе 4 121
Заключение 122
Список литературы 124
- Анализ существующих методов восстановления изображений земной поверхности, полученных в оптическом диапазоне
- Адаптивное оценивание параметров восстанавливающего оператора
- Исследование алгоритма восстановления на модельных примерах
- Проведение сравнительных вычислительных экспериментов
Анализ существующих методов восстановления изображений земной поверхности, полученных в оптическом диапазоне
Пространственная разрешающая способность большинства современных космических систем формирования изображений определяется количеством фотоприёмных элементов ПЗС-матриц. Улучшение этой характеристики путем изготовления ПЗС-матриц с большим числом фотоприёмных элементов и создания на их основе новых оптоэлектронных приборов является наиболее очевидным, и, в тоже время, дорогостоящим и неоптимальным методом [100].
В последние годы в технике получения высококачественных космических изображений, в том числе с поддержкой высокого и сверхвысокого разрешения, наметилась тенденция широкого применения методов цифровой обработки пикселей, поступающих из фотоприёмных элементов [33, 34]. Получаемые в этом случае после обработки цифровые пиксели изображений – т.н. субпиксели – имеют гораздо меньшую апертуру, чем пиксели первичных данных без цифровой обработки. Физический смысл таких методов заключается в том, что увеличение разрешения и снижение апертуры пикселей на одном и том же кадре изображения связано с необходимостью получения реальных дополнительных пикселей изображения исследуемой поверхности с последующим уменьшением их размера и размещением их с высокой степенью точности на результирующем транспаранте по методу пересечения взаимно смещенных в пространстве одинаковых решетчатых функций [1, 101]. Изображение объекта попадает на соседние ПЗС-линейки, создающие эффект избыточной дискретизации, на основе которой решением системы уравнений находится новая система пикселей, поддерживающая улучшенное почти вдвое разрешение, по сравнению с формированием изображения на одиночных линейках ПЗС. При этом применяется не тривиальная интерполяционная процедура вычисления промежуточных пикселей, а реализуется увеличение числа точек, дискретизирующих исходное изображение и, в соответствии с теоремой Котельникова [35], создающее более широкую пространственно-частотную передаточную характеристику системы формирования изображения [36,37,38] -частотно контрастную характеристику (ЧКХ).
При этом в процессах улучшения ПР, особенно в этапных или итеративных, не всегда достаточно точно можно указать момент необходимого перехода к субпиксельному представлению изображения. В общем виде снижение апертур ФРТ на изображении и подавление шумовой составляющей проводится на реальном изображении или его модели со спектром пространственных частот, пересекающимся с полосой аналогичного спектрального портрета шумовой составляющей изображения. Такая модель является актуальной, т.к. радиусы корреляции контуров, вычисляемые поперек простирания контуров объектов и радиусы корреляции шумов совпадают по значению в силу ограничения пространственно-спектральной полосы всего тракта дистанционного зондирования [39].
В условиях, когда функция рассеивания точки имеет практически нулевой диаметр, формула (1.1.1) дает возможность определить размеры пикселя на Земле при заданных размерах пикселей на ПЗС-матрице фотоприёмника или поставить задачу о достижении необходимого размера пикселей на ПЗС-матрице при требуемых размерах пикселей на наблюдаемой поверхности [21, 22, 39, 40, 41, 42, 43, 44, 45, 102].
Отталкиваясь от определения разрешения по Рэлею [47] с линейной метрикой (или выполняя обращение метрики по Фуко [25]) с учетом формулы (1.1.1) априорное ПР (ПРа) можно представить как отношение энергии сигнала, соответствующей в данной разрядной сетке представления изображения скачку кода на единицу в младшем разряде к энергии сигнала максимального значения кода цифрового представления изображения.
На практике (например, при получении изображений спутником Ресурс-Д) используется апостериорное ПР (ПРап), измеряемое отношением энергии полезного сигнала к энергии шума на изображении после подавления шумов квантования, атмосферных, аппаратных искажений, а также шумов, связанных с несовершенством процедур коррекции и предварительной обработки формируемого изображения [22]. Известно, что различные процедуры восстановления изображения подавляют шумовую и высокочастотную составляющие изображения. Это подавление можно с успехом использовать для подготовки полезной и шумовой составляющих при оценках ПРап. При этом сохранение структуры изображения, т.е. границ структуры, представляется весьма ценным [48].
Несмотря на некоторую неоднозначность в определениях и неустановившуюся еще терминологию в определении понятия «сверхразрешение», которым нередко пользуются для определения методов формирования изображения с меньшими расстояниями между его элементами (по Рэлею), чем позволяет современная промышленная технология преобразования данных, целесообразно в данной работе считать сверхразрешением перепикселизацию с уменьшением апертуры пикселя в соответствии с теоремой Котельникова.
В работах по сверхразрешению на формируемых изображениях исследованы альтернативные критерии оценки четкости и дешифрируемости объектов на изображении, основанные на расчете энергии высокочастотных составляющих спектра изображения, дисперсии разности соседних элементов и расчете энтропии. Экспериментально установлено, что эти критерии дают близкие относительные оценки качества, а при комплексировании данных от двух субпиксельно смещенных ПЗС-матриц можно достичь повышения разрешающей способности съемки в 1,5 раза, а при комплексировании данных от трех ПЗС-линеек – в 1,9 раза.
В целом данное направление работ можно отнести к технике коррекции изображений и поддерживаемой разрешающей способности с помощью улучшения и превышения штатных возможностей аппаратуры [44, 49, 103].
Применение субпиксельных представлений оказывает сильное влияние на организацию процедур коррекции изображений. Например, при смазе изображения за счет ухода скорости движения изображения (СДИ) от номинала часто возникает нецелое отношение апертур смаза к расстоянию между строками формируемого изображения, что при использовании стандартных пиксельных процедур коррекции реализует достаточно приближенный результат.
Адаптивное оценивание параметров восстанавливающего оператора
В данном разделе представлена разработка синтеза частотно-контрастной характеристики тракта формирования для сформированного космического изображения. В разделе 2.1 построено представление повышения разрешающей способности с использованием оператора (l+agrada) с дифференциальными операторами нецелого порядка. Положительная определенность нормы и операторов, определенных в пиксельном пространстве изображений, диктуется необходимостью инверсии результирующих и присваиваемых пикселям изображений дополнительных кодов, получающихся на неположительных значениях вычисляемых яркостей в монохромных каналах.
Оператор (l+agrada\ описанный в разделе 2.1, не является средством прямого решения уравнения Фредгольма - уравнения формирования наблюдаемого изображения по восстанавливаемому изображению, а служит для синтеза ЧКХ, в последствии оптимизируемой. Такую ЧКХ можно получить, используя как само изображение, так и спектральное представление (СП) оцененной на изображении ФРТ. Это СП как правило используется в степени -1 как СП ядра оператора деконволюции, модернизируемого регуляризацией и дополнительными членами по имеющейся априорной информации, например, в формуле Винера-Тихонова [42, 64].
Обширная теория и практика высококачественного восстановления и фильтрации изображений оперирует с оптимальными ЧКХ с квазипрямоугольной огибающей на любом азимутальном срезе, т.е. с классическим окном пропускания пространственных мод СП. Подбор параметра СП имеющегося оператора, согласованные с оптимальной ЧКХ тракта, порождают решения по восстановлению и могут быть выбраны пользователем из соображений выдерживаемой погрешности процессов, целевых уклонов при оптимизации ЧКХ и дополнительных требований к оптимальности [25, 65]. Для тракта ДЗЗ необходимо знать значение верхней моды сов в его ПЧС. В настоящей работе операторы А конволюции и деконволюции -интегральные; ограниченные; определены для элементов х линейного метрического Гильбертового пространства с нормой [66]:
Изображения дискретизируются в соответствии с теоремой Котельникова, являются элементами х конечномерного линейного, метрического Гильбертового пространства с евклидовой нормой; СП интегральных операторов (интегралов свертки) задаются проекциями их ядер на ортонормированный базис тригонометрических функций. В пространстве, сопряженном с , при помощи дискретного преобразования Фурье определяется градиентная операция как норма (порождающая на разности векторов метрику) на дифференциальных операторах вида (2.1.9).
Оптимизация аддитивного представления [62, 65], также как и ЧКХ, заключается в нахождении значений параметров а и а в операторе (1+a(grada)) с получением сопутствующей ЧКХ, удовлетворяющей набору объективных требований с учетом СП вычисляемой функции смаза от опознаваемой на космическом изображении границы (функции рассеяния границы (ФРГ)) или линии (функции рассеяния линии (ФРЛ)) [1] или фрагмента полигонного ориентира. Это реализуется на космическом изображении, как отмечалось выше, без учета кинематических особенностей орбитальной съемки, которые меняют результаты коррекций по расчетам и измерениям ФРТ [16]. Для любого азимутального среза оптимальная ЧКХ реализуется по огибающей как квазипрямоугольное окно пропускания мод спектральных представлений входного изображения. Сохранение характерного для трактов ДЗЗ спада ЧКХ приводит к недостаточнойкомпенсации действия ФРТ, а превышение его уровня компенсации - к усилению контраста и стиранию сверхмелких деталей изображения при нормализации результатов деконволюции [104].
Отметим свойства гипотетической оптимальной ЧКХ при восстановлении: на всех этапах деконволюции все ЧКХ совпадают в низко- и среднечастотной части, и только в высокочастотной части различаются; все ЧКХ, соответствующие недостаточному восстановлению вписываются в оптимальную ЧКХ; ЧКХ с подчеркиванием высших мод или с контрастированием и, тем более, выделением контуров мажорируют все ЧКХ; оптимальная ЧКХ не противоречит спектрам ФРТ, ФРЛ или ФРГ от полигонного артефакта [62].
Исходя из вышесказанного, предлагаются следующие условия условия: а) максимизации функционала - объема под огибающей ЧКХ (соотношение 2.2.3); б) отсутствие тенденции спада и роста аппликат ЧКХ вплоть до приближения к высшим спектральным модам тракта - примерной гладкости и параллельности огибающей поверхности ЧКХ координатной плоскости (соотношение 2.2.4, 2.2.5,соответственно). Это условие имеет всегда максимальный приоритет и, соответственно, высокое значение весового коэффициента при вычислении оптимального а и а ; в) обеспечение начала устойчивого спада огибающей ЧКХ с некоторой моды со2 сов, причем с минимизацией разницы сов - со2 (соотношения 2.2.6, 2.2.6 а, 2.2.7, 2.2.8); г) вписывание огибающей спектрального представления, наблюдаемого на изображении линеамента или границы в огибающую ЧКХ, включая вписывание в огибающую ЧКХ СП, обратного результату применения фильтра (соотношение 2.2.9); д) непревышение значением а порога 0,25, при котором присутствуют глобальное контрастирование или выделение контуров [65].
Следует отметить, что одно из самых естественных требований к ЧКХ, а именно - инвариантность ЧКХ к повороту, может относиться только к ЧКХ на выходе оптоэлектронного преобразователя. Инвариантность нарушается для всего тракта формирования изображения, т.к. в тракт включаются сканирование строк чувствительных матриц и внесение коррекций на угол поворота оптической оси системы, на тонгаж и рыскание аппарата.
Исследование алгоритма восстановления на модельных примерах
Лабораторные математические модели представления изображений с улучшенным разрешением на основе субпиксельных сдвигов и ряда предварительных расфокусировок, проецируемых на приборы с зарядовой связью (ПЗС) изображений, показывают неоправданно обнадеживающие результаты, тогда как реальные достижения измеряются показателями улучшения в 1,4-1,8 раза [71]. Важным фактором при этом является геометрическое несовершенство элементов и строк пикселей ПЗС, а так же соотношение между радиометрическим разрешением и средним порогом маскировки шумом огибающей, сформированной при расфокусировке ФРТ. Что касается комплексирования улучшенного изображения из фрагментов ранее отснятых снимков, то выполнение требований по нивелированию невязок фрагментов или практически невозможно или требует неоправданно высоких вычислительных затрат [25,102].
При этом совмещении одного из выбранных фрагментов к состоянию необходимого сдвига на долю пикселя, требуемого субпиксельной технологией, реализуется интерполяцией промежуточных отсчетов на изображении с последующей децимацией и сохранением отсчетов, сдвинутых относительно исходных. При выполнении интерполяции методом Уиттекера-Шеннона осуществляется достаточное снижение среднего значения распределения вероятности паразитного сглаживания фронтов градаций яркости на изображениях (в том числе и при выполнении промежуточных вычислений), выполняемое интерполяциями по соседним отсчетам и другими видами интерполяции, что негативно сказывается на расширении радиуса пространственного спектра результирующего изображения получаемого субпиксельными технологиями. При этом используется свойство сохранения ширины спектра интерполируемой функции методом Уиттекера-Шеннона: возможность расположить вычисляемый дополнительный отсчет достаточно близко к одному из исходных отсчетов, чем при итоговом формировании интерполированного изображения. Это позволяет сохранить (не сгладить) имеющийся заметный скачек градации яркости (в общем случае критерии наличия скачков - в исходной пиксельной системе вероятностные и могут подключаться к процедуре интерполяции в автоматизированном режиме анализа исходного изображения). Соблюдение этих правил при использовании, например, интерполяции по расположенным рядом пикселям с учетом пропорционального деления интервала не сохраняет ширину спектра исходного изображения [74].
На основе моделирования влияющих негативно факторов предложена схема формирования субпиксельного изображения, дающая эффективность улучшения в пределах 20%, организованную на трех ПЗС-линейках (или матрицах) с совмещением в схеме и приемлемого уровня расфокусировки за счет различия апертур пикселей от одного ПЗС к другому. Техника субпиксельного сдвига реализуется на сдвинутых относительно друг друга идентичных линейках на треть апертуры пикселя, но с подачей на вычислитель полусуммы каждой соседней пары элементов одной ПЗС и трети суммы с каждых трех соседних элементов другой ПЗС-линейки. Для двух субпиксельно сдвинутых линеек эффект улучшения разрешения, как установлено эмпирически, несколько превышает значение 4Ї, для трех линеек - V3 . Однако, перспективность такой схемы субпиксельных технологий (сверхразрешения) на практике невысокая, т.к. нагромождение ПЗС не двух, а уже трех (на каждый из спектральных каналов) потребует ужесточения к требованиям по реализации оптической системы, увеличения габаритов БОЭА, дополнительных средств и т.п. [73].
Стандартная (редуцированная) схема субпиксельного сдвига для формирующих сверхразрешение строк (ПЗС или изображений) Алгоритм P(xi) формирования системы уравнений для нахождения субпикселов xi, соответствующих выполненной таким образом передискретизации и восстановления потерянной на исходных строках тонкой структуры
Система уравнений переопределена и позволяет составить 13 уравнений. Построим математическую модель влияния фактора геометрического несовершенства ПЗС на снижение радиуса пространственно-спектрального портрета изображения, восстанавливаемого в субпикселях на двух взаимно сдвинутых ПЗС.
Очевидно, что сформированная случайными поперечными сдвигами в позициях пикселей ПЗС ФРТ при режиме в силу теоремы о линейном преобразовании случайной величины преобразуется в ФРТ с гауссовой огибающей [38,69,90]. При организации технологии сверхразрешения кроме ФРТ т.н. стандартного типа (превышающих пиксели по апертуре) к формированию изображения присоединяется и сверхмалая ФРТ, поэтому целесообразно реализовать модель на поперечных взаимных сдвигах [74].
Это соответствует наличию на изображении, полученном в режиме сверхразрешения неучтенной, а главное трудно измеряемой ФРТ, которая занижает КПД субпиксельной технологии, т.к. показатель степени при экспоненте в выражении (3.4.5) ненулевой (в силу технологических причин), хотя достаточно мал (дисперсия на порядок меньше апертур самих пикселей).
Эти соотношения связывают разрешение по координатам с разрешением по частотам, а также разрешение на изображении с количеством отсчетов спектра изображения. Каждый отсчет дискретного спектра представлен двумя составляющими (ReFs и ImF5), но спектр вещественной функции S(r) удовлетворяет соотношению Fs(-a ) = F(a ), т.е. и функция, и ее спектр содержат равное число независимых отсчетов и информационно равноценны: Fs{-сох-соу) = F(cox,coy); здесь F - комплексно-сопряженная величина.
Выделяемая двумерная F(PSF) в двумерной модели может не быть гауссовой поверхностью и сформированный тракт не позволяет решать задачу повышения разрешения на SR, так как спектральный портрет SR строго ограничен модой сов. Таких огибающих ЧКХ может быть достаточно много: полуволна синусоиды, колоколообразная функция, прямоугольная, трапецевидная и т.д. Главное здесь то, что гауссоподобная ЧКХ никогда не пересекается с осью частот, что определяет, в принципе, значительную вероятность постановки и решения задачи восстановления. В этом случае подавленные пространственные гармоники восстанавливаемого изображения 57 можно усилить подбором операций дифференциального характера (градиентных фильтров с порядком а со значением не выше 0,25). Спектральные портреты, которых являются неубывающими функциями от (Jco)a\. С физической точки зрения это означает, что можно применить процедуру восстановления, описанную в диссертационной работе.
Известно, что ПЧС каждой строки паттерна, предназначенного для формирования сверхразрешения, имеет вид в каждой строке Fp(aj) = F(aj)exp(-aj2a2r) (3.4.9) В формуле (3.4.9)учитываются составляющие от алгоритмов бортовой коррекции при формировании изображения. Значение а г, при условии всех выполненных коррекций индустриального типа, определяет с высокой степенью вероятности значением сов, Предположим, что и вдоль столбцов справедлива формула (3.4.9).Так как угловые повороты орбитального аппарата, крен и тонгаж, а главное, изменение разности реальной и программируемой высоты при полете из-за неучтенной или недостаточно детальной модели рельефа местности, то вносимые изменения в методе, выполняемые в этих случаях линейными преобразованиями, т.е. весовыми суммированиями с нормализацией, практически аналогичны описанной выше модели. Для ПЧС входных в технологию сверхразрешения четырех изображений запишем соотношение
Проведение сравнительных вычислительных экспериментов
Решение уравнения Фредгольма, соответствующего модели Бейтса и Мак Доннела, а для формируемого космического изображения относительно & это деконволюция, относится к некорректно поставленным задачам, т.к. малые вариации в SR могут приводить к большим изменениям в Si. При этом используются две группы методов деконволюции - прямые и итерационные. В прямых методах для стабилизации решения используются различные способы регуляризации, самым распространенным из которых является регуляризация Тихонова, определяющая решение Si л, которое минимизирует выражение: SiX = argmin (Si){\\F1(H) SrSR\\2 l\\L(Si) \\2} (4.2.1) где первое слагаемое - невязка, L - соответствующая норма решения, а параметр регуляризации X задает вес минимизации нормы решения по отношению к минимизации нормы невязки. Большие X отбирают решения с малой нормой, но с большей невязкой. Малые X соответствуют решениям с большой нормой, но малой невязкой (в итерационных методах деконволюции роль параметра регуляризации выполняет номер итерации). Итерационный процесс можно наблюдать в пространстве решений, и, таким образом, эти методы позволяют избежать трудностей, связанных с выбором оптимального параметра регуляризации. При этом множество решений, даваемых любым методом деконволюции, требует критерия отбора оптимального решения [74].
Отметим соотношение между методом регуляризации Тихонова и методом восстановления с одновременным повышением резкости космических изображений высокого разрешения, требованием к сохранению валидности которых является именно восстановление с повышением резкости, т.е. компенсация сглаживающего действия ФРТ, а не привнесение в изображение чрезмерного усиления высших мод его пространственного спектра.
В предложенном методе аддитивное представление восстанавливающего оператора с оптимизацией сопутствующей коррекции ЧКХ тракта формирования изображения строится при условии отсутствия провоцирования выделения контуров или процессов контрастирования, зависящих от вида H, определяемого на опорном ориентире. Реализуется ближайшее по заданной погрешности приближение ЧКХ к прямоугольному виду на любом азимутальном срезе. Единственным источником множественности решения при этом является выбор весовых коэффициентов для оптимизации ЧКХ.
Приведенное выше описание регуляризации Тихонова при любом выборе соотношений между невязкой и регуляризационным параметром несвободно от скрытой в СП ядра деконволюции возможности впечатывания контуров в изображение или чрезмерного повышения высших мод ЧКХ. Контроль по ЧКХ в методе не предусматривается и схема метода этого не позволяет. Единственным возможным вариантом сохранения измерительных свойств изображения является тщательный выбор опорного ориентира со спектральным представлением ФРТ таким, чтобы не появлялись резидентные явления, что с технической точки зрения не всегда реализуемо.
Это же можно сказать и о таких популярных и достаточно весомых при обработке космических изображений методах итеративной деконволюции без привлечения баллистических параметров и смежных данных из каналов тракта ДЗЗ как метод Ван Циттера (на нем в диссертационной работе построена коррекция резкости с возмущением оператора деконволюции), метод Джанссона, метод Голда, метод максимума правдоподобия (ММП), метод разложения по сингулярным значениям [82,83,84].
Отметим, что ММП практически повторяет формулу метода Джанссона, но в нормированном варианте изображений и ФРТ с определением шагового параметра максимизацией функции правдоподобия, т.е. минимизации разницы на итерациях между искомым решением и решением в предыдущих итерациях.
В практике выполнения деконволюции описываемые методы ранжируют по степени повышения качества решения в следующем порядке: 1) метод Ван Циттерта; 2) метод Джанссона; 3) метод Тихонова; 4) метод Голда; 5) ММП. Однако при этом в результате работы методов превентивные меры по компенсации усиления контраста не предъявляются, что делает соответствующие оценки качества улучшения резкости недостаточно корректными.
Ниже приведены результаты восстановления с повышением резкости изображения высокого разрешения с аппарата OrbView-3 методами Ван Циттера (без оптимизации ЧКХ), Джанссона, Голда и метода восстановления, разработанного в данной работе. Опорный ориентир выбирался на модели изображения, полученной методом расфокусировки на стандартном графическом пакете [75].
Следует отметить, что соотношение (4.2.1) показывает квазилинейное нарастание значения минимума в зависимости от уровня восстановления изображения. Такое возрастание значения минимума выражения (4.2.1) проявляется при обработке изображения, когда оно преобразуется к изображению с резко подчеркнутыми контурами и когда подавляется низкочастотная часть пространственно-частотного спектра. При этом (в эксперименте) реализуется итеративный процесс вида (2.3.6): F(0)(SI)=F(SR), . . F(n)(SI)=F(SR)+(1-H)F(n-1)(SI), остановка в итеративном процессе осуществляется оператором – интерпретатором. Каждый F(n)(SI) отличается от F(0)(SI) увеличением амплитуд верхних мод ПЧС с увеличением номера n. Единственным средством избежать усиления контраста на изображении в процессе его восстановления является задание ограничения по изменению уровня наклона возрастания верхних мод ПЧС при обработке относительно исходного изображения. Это эквивалентно контролю итераций по сопутствующей ЧКХ или контролю по уровню увеличения эффективного порядка соответствующего оператора при аддитивном представлении.