Численное решение уравнений в свертках и алгоритмы деконволюции
Искажение сигналов и размытие изображений
В физике, астрономии, оптике и цифровой обработке сигналов мы постоянно сталкиваемся с тем, что измерительная аппаратура несовершенна. Любой прибор (будь то телескоп, микроскоп или сейсмограф) неизбежно искажает истинный физический сигнал. Математически процесс искажения линейной стационарной системой (LTI-системой) описывается операцией свертки: истинный сигнал «сворачивается» с так называемой аппаратной функцией прибора (функцией рассеяния точки, PSF). В результате мы получаем размытую, сглаженную картину реальности.
Обратная математическая задача — восстановить истинный, первоначальный сигнал по искаженному результату измерений и известной аппаратной функции — называется задачей деконволюции (обратной свертки). С математической точки зрения, деконволюция сводится к решению интегрального уравнения Фредгольма первого рода со специфическим ядром, зависящим от разности аргументов. Эта задача является классическим примером некорректно поставленной проблемы: малейший шум в измерительных данных (а он есть всегда) при попытке прямой деконволюции приводит к катастрофическому разрушению результата, топя истинный сигнал в высокочастотных осцилляциях.
Деление в частотной области и фильтр Винера
Благодаря теореме о свертке, сложная интегральная операция свертки в физическом времени (или пространстве) превращается в простое алгебраическое умножение в частотной области (спектральном пространстве). Соответственно, кажется, что деконволюцию можно выполнить элементарно: применить Быстрое Преобразование Фурье (БПФ) к измеренному сигналу и аппаратной функции, поделить первый спектр на второй, а затем применить обратное БПФ. Этот метод называется обратной фильтрацией.
Однако на практике метод прямого деления спектров абсолютно не работает. Высокие частоты аппаратной функции обычно стремятся к нулю (прибор срезает резкие скачки). При делении зашумленного спектра на числа, близкие к нулю, мы получаем гигантские значения (умножение шума на бесконечность). Для стабилизации алгоритма Норберт Винер разработал оптимальный стохастический фильтр (фильтр Винера). В спектр знаменателя математически добавляется специальное слагаемое, пропорциональное отношению мощности шума к мощности полезного сигнала. Этот регуляризующий член не дает знаменателю обратиться в ноль, подавляя частоты, где доминирует шум, и пропуская частоты, где доминирует полезный сигнал. Фильтр Винера до сих пор является золотым стандартом для восстановления расфокусированных цифровых фотографий и очистки старых аудиозаписей.
Итерационная деконволюция Ричардсона-Люси
В астрономии и медицинской микроскопии, где измеряются интенсивности света (которые по физическому смыслу не могут быть отрицательными), фильтр Винера дает сбои, так как часто порождает отрицательные «волны» (артефакты звона) вокруг ярких объектов. В 1970-х годах Уильям Ричардсон и Леон Люси независимо друг от друга вывели мощный итерационный алгоритм деконволюции, основанный на теореме Байеса и максимизации функции правдоподобия в предположении пуассоновской статистики фотонного шума.
Алгоритм Ричардсона-Люси работает итеративно: на каждом шаге он берет текущее приближение истинного изображения, сворачивает его с аппаратной функцией (симулируя искажение), сравнивает полученный результат с реальной размытой фотографией путем деления, а затем использует эту разность для коррекции приближения на следующем шаге. Гениальность метода состоит в том, что он математически гарантирует неотрицательность пикселей на любой итерации и отлично сохраняет полную энергию сигнала. Именно алгоритм Ричардсона-Люси стал спасением для космического телескопа «Хаббл» в 1990 году, когда выяснилось, что его главное зеркало отполировано с ошибкой (имеет сферическую аберрацию): численные методы деконволюции позволили ученым получать резкие снимки галактик еще до прибытия ремонтной экспедиции с корректирующей оптикой.
