Бойков И.В., Мойко Н.В. ОБ ОДНОМ МЕТОДЕ ВОССТАНОВЛЕНИЯ ИМПУЛЬСНОЙ ФУНКЦИИ И ИСКАЖЕННОГО АТМОСФЕРОЙ КОРОТКОЭКСПОЗИЦИОННОГО ИЗОБРАЖЕНИЯ
Введение.
Изображение объекта, наблюдаемого через турбулентную атмосферу, является сильно искаженным. Для восстановления изображения существуют различные методы фильтрации [1-3]. В [1,2] излагаются методы фильтрации, требующие знания оптической передаточности функции (ОПФ) Н (ш) или импульсной функции
Н(£) . Так как априори такая информация, как правило, отсутствует, то возникает задача одновременного
определения импульсной функции и входного сигнала.
В [3] изложен метод решения указанной задачи, основанный на осреднении фазы ОПФ. Ниже предлагается метод, основанный на теории краевых задач Римана [4,5].
Следуя [3], задача формулируется следующим образом.
Имеется короткоэкспозиционное изображение I м протяженного объекта
I (м ) = /х (м) Н(м - м ) (0.1)
где х (м)—неизвестное истинное распределение интенсивности объекта со сложной конфигурацией, а Н (м) случайная импульсная функция системы "атмосфера-телескоп".
Предполагается известным некоторый функционал от функции Н(м). В качестве такого функционала мож-
но взять
да
| Н (м) ds = к
—да
или
да
| |Н (м)| ds = к,
—да
где к — еот1,1 = 1,2,....
Можно также считать известной функцию Н(м), полученную по некоторому алгоритму осреднения функции н (м).
Предполагается, что в некоторой области, размеры которой меньше чем размеры области изображения, Н (м ) > 0, а в остальных точках области изображения Н (м ) = 0.
Уравнение (0.1) имеет бесконечное число решений и поэтому для одновременного восстановления входного сигнала х (м) и импульсной функции Н (м) требуется дополнительная информация. В качестве такой
информации можно задать значения функционалов от х (м), н (м) или от их преобразований Фурье.
Ниже будем считать заданным значение функционала
J | h (s) | ds = к,
(0.2)
где k =const - постоянная.
Как будет видно из дальнейшего, вид функционала не влияет на проводимые выкладки, и полученные результаты справедливы для любых функционалов.
Ставится задача: найти единственные значения функций x (s) и h (s), удовлетворяющие уравнению (0.1)
и условию (0.2).
Поставленная задача может быть решена двумя способами.
По первому способу изображение f(x) разбивается на фрагменты ширины d, где d полагается достаточно малым; задача восстановления изображения на каждом фрагменте сводится к решению одномерного уравнения в свертках (см. ниже (1.3)). Метод решения этого уравнения изложен в параграфе 1.
В результате получаем набор решений x. (s), h (s), j = 1,2,..., w, где n - число выделенных фрагментов.
Затем из изображения функции f (s) выделяется фрагмент шириной d со сторонами, параллельными другой координатной оси, составляется и решается одномерное уравнение. Решение этого уравнения
(s(s2 ), h(s2 )) является «огибающей», определяющей коэффициенты, на которые умножаются значения xj (si),hj (si), j = 1,2,..., w-
В результате этих действий восстанавливаются функции x(s^,s2 ), h(s1 , s2). Единственное решение
определяется из условия (0.2).
По второму способу входной сигнал x (s) восстанавливается как решение краевой задачи Римана в би-
цилиндрической области. Этот метод изложен в параграфе 2.
§ 1. Первый способ сведения к одномерным уравнениям.
1. Вначале исследуем теоретическую возможность одновременного восстановления аппаратной функции h(s) и входного сигнала x(t) по известному выходному сигналу f (t) , равному
| Н (г — г) х (г) dт = I (г), — да< г <да. (1.1)
—а
Предполагается, что функции Н (г) и х (г) финитны с носителем [— а,а] и суммируемы с квадратом. Отметим, что симметричный относительно начала координат сегмент [—а,а] взят только из удобства обозначений.
Поэтому уравнение (1.1) можно представить в виде
да
| Н(г — г)х(г)dг = I(г), —да<г<да. (1.2)
—да
Сделаем в уравнении (1.2) следующие замены переменных: г = V — а, г = ш + а. Т о г д а уравнение (1.2) име-
ет вид
да
| Н(V — ш — 2а)х(ш+а)dш = I(V — а), —да<V<да (1.3)
—да
Введем обозначения: Н (г )= Н (г — 2а), х (г ) = х (г + а), I (г )= I (г — а).
Запишем уравнение (1.3) в этих обозначениях:
0
| Н1 (V — ш)х1 (ш)dш = Ц(V), —да< V <да. (1.4)
—2а
Обозначим через Н (ш), Х^ (ш), (ш) преобразование Фурье функций Ну (г), х^ (г), I (г).
Применим к уравнению (1.4) преобразование Фурье. Так как функция кх (г) равна нулю на полуоси (—да,0), а х (г) - на полуоси (0, да) и обе эти функции финитны и суммируемы с квадратом, то ([5, стр.
25]) функция Н(ш) будет аналитической в верхней полуплоскости, а функция Х1 (ш) - в нижней полу-
плоскости.
После применения преобразования Фурье, уравнение (1.4) принимает вид Н+ (ш)Х1“(ш)= ^ (ш). (1.5)
Так как ^ (да) = ^ (—^) = 0, то уравнение (1.5) относится к исключительному случаю краевых задач.
Пусть функция ^(ш) имеет нуль порядка г на бесконечности.
Тогда уравнение (1.5) можно представить в виде
Н+ (ш)Х[ (ш) = (ш + I)г F2 (ш)
или, введя новые аналитические функции
О+ (ш) = Н+ (ш) / (ш + 1)г, О- (ш) = Х~ (ш),
к краевой задаче
О+(ш)О (ш) = F2 (ш). (1.6)
Таким образом, задача свелась к решению уравнения (1.6), которое может быть получено факторизацией функции F(ш). Факторизация функции F(ш) несколько отличается от стандартной факторизации [4,5], но тем не менее уравнение (1.6) сводится к видоизмененной задаче о скачке.
Пусть индекс х функции F (ш) равен 0. В этом случае 1пF2(ш) будет однозначной функцией. Логарифмируя уравнение (1.6), имеем 1пО+ (ш) + 1п О~ (ш) = 1пF2 (ш).
Введем функцию
у (г) = — | 1пF (ш) е ‘шгdш. (1.7)
2ж
Тогда решение краевой задачи (1.7) представляется в виде [5]
о± (г )=ег±(г),
где
1 да
Г+( 2 ) = -^= |у(г )sgn e,zГdГ, (1.8)
0 0
(1.9)
—
Таким образом, при х=0 получено единственное (с точностью до постоянного множителя а : Оа+ (2) = аО+ (2 (, С+ (2 ) = а — (2 ) ) решение задачи одновременного восстановления аппаратной функции и входного сигнала.
В случае, когда х>0, краевая задача (1.6) имеет х линейно независимых решений [4,5].
В случае, когда х<0 задача (1.6) не имеет решения [4,5].
Приступим теперь к изложению алгоритма одновременного восстановления входного сигнала и импульсной функции по уравнению (0.1) и дополнительному условию (0.2).
Будем для определенности считать, что входной сигнал / (м, м2) имеет носитель [а, Ь;с, d ]. Введем узлы Ук = с + (^ — с)к/п, к = 0,1,...,п. Пусть Н = (Л — с)/п.
Каждому узлу ^ поставим в соответствие уравнение Ь
|Нк (м — а) хк (а)^= /к (м^
а
где Л (м) = /(s,v^),k = 0,1,...,п — 1.
Решив эти уравнения, получаем П наборов решений (х^ (м),Н (м)), к = 0,1,...,п — 1.
Зафиксируем теперь значение ш = (а + Ь)/2 и рассмотрим уравнение d _ _
1Н (м — а) х(а)dа = I (ш,м).
с
Решив эти уравнения, выберем одно из решений, скажем (н (м),х (м)).
Теперь согласуем наборы решений (хк (м),Ик (м)), к = 0,1,...,п — 1 с решением (х (м),Н (м)). Для этого, полагая х (уу) * 0,у = 0,1,...,п — 1, введем в рассмотрение новый набор решений (хк (м),Нк (м)), где
хк (м) = хк (м) х (vk )/ хк (^) при хк (vk )* 0, хк(м )=0 при хк (vk )=o,
нк (м)=нк (м) ^)/Нк (vk) при к ы* 0 к=o,l,..., п—1,
\ (м) = 0 при Ик (м) = 0, к = 0,1,..., п — 1.
Замечание 1. Если при каком - нибудь значении V хк (^) = 0, то вместо того, чтобы полагать хк (м) = 0, естественно считать х^ (м )=( хк+1 (3)+хк—1 (м))/2.
Решение уравнения (0.1) определяется как функции (х(м, м2 )) где
х(м,^) = Ахк (м),Н(м,м2) — А 1Ьк (м) при уk < м2 < Уk+1,к = 0,1,...,п — 1. Константа A выбирается из условия (0.2).
Замечание. При окончательном определении функций х(м, м2) и Н(м, м2) естественно провести их осреднения воспользовавшись ядрами усреднения с диаметром усреднения равным 3Н,4Н, 5Н.
§2. Второй метод. Аналитическое решение многомерных уравнений.
Исследуем возможность восстановления изображения конечного плоского объекта произвольной достаточно сложной конфигурации. Задача восстановления изображения в данном случае описывается уравнением
а а
| |к(гх — гг2 — г2)х(г,г2^г^г2 = I(гх,г2). (2.1)
—а —а
Будем считать, что функции Н(г1,г2 X х (А г) имеют носитель [—а, а]2, функция I (^ ,г2) определена в квадрате [—2а, 2а]2.
Так как функция х (г ,г2) финитна в квадрате [—а, а]2, а функция I (^, г2) может быть доопределена нулем вне квадрата [—2а,2а]2, то уравнение (2.1) можно представить в виде
да да
| | Н (г1 —г1г2 —г2) х (г1,г2) dг1dгг = I (г1, г2), — да< г1, г2 <да (2.2)
—да —да
Введем новые переменные V = Г+ а,V =Г+ +,ш = Ч — 2а, ш2 = г2 — 2а.
В результате подстановки предыдущее уравнение (2.2) принимает вид
да да
| | к(ш1 — V + 3+,ш2 — v2 + За)х(V — а,v2 — а)dv1dv2 = I(шг + 2а,ш2 + 2а). (2.3)
—да —да
Введем функции
Н(ш,ш) = Н(ш + За,ш2 + За), х(V,V) = х(V — а^г — а), I(ш,ш) = I(ш + 2а,ш2 + 2а).
Уравнение (2.3) преобразуется к виду
1 1 Н (ш>1 — х,ш2 — v2 ) х1 (vl, v2 ) dVldV2 = I (vl, V ). (2.4)
СО —СО
Обозначим через гк, к = 1,2, плоскость комплексной переменной = ш+ V.
Через В+(В—) обозначается верхняя и нижняя полуплоскость плоскости 2к (к = 1,2,) . Символы
В++, В+ , В +, В обозначают топологическое произведение областей
В±(к = 1,2):В++ = О+ В+, В+— = В+В—,В_+ = В“В+,О" = В— В~. Через ^++(ш,ш), ^+— (ш,ш2), ^_+(ш,ш), ^“(шш)
обозначаются предельные значения функций ^++( 2, 22 ),^+ (2,2 ),^ +(21,22 ),^ (2,22 ) аналитических соответ-
ственно в областях В++, В+ , В +, В при стремлении точек (, 22) к остову шх®2 по любым направлениям не касательным к осям ш или ш2. Так как по предположению функции Н (гх ,г2) и х (гх ,г2) финитны с носителями расположенными в квадрате
[—а,а]2, то функция К (Ч,г2) отлична от нуля только в первом, а х1 (г1, г2 ) только в третьем квадрантах. Следовательно, применение преобразования Фурье к уравнению (2.4) приводит к краевой задаче Римана Н++ (ш, ш2 )X— — (ш, ш2 )= F (ш, ш2 ). (2.5)
Пусть функция F (ш,ш) представима в виде
F (шш) = (ш —О*1 (ш —0^ ^ (ш,ш),
причем функция ^ (ш,ш) не обращается в нуль на остове ^х^ и ее частные индексы [6] равны нулю.
Тогда краевая задача (2.5) эквивалентна
О++(ш,ш)О (ш,ш2) = F2 (ш,ш), (2.6)
где
О++ (ш, ш) = Н^ + (ш, ш2), о~ _(ш,ш) = (ш —О*1 (ш —0^ х^ — (ш,ш).
Аналитические методы решения уравнений вида (2.6) изучались в [6]. Следуя [6] укажем один из способов решения уравнения (2.6). Так как функция F2(ш,ш2) не обращается в нуль на остове ш^хт2 , и ее
частные индексы равны нулю, то существует функция Ф(ш,ш2) = 1пF2 (ш,ш2). Применим к функции Ф(ш,ш) формулу Сохоцкого [4,стр.91],
ф(ш,ш)=ф++ (шш)—ф+— (ш,ш)—ф—+(ш,ш)+ф—— (ш,ш),
где
ч____^ да Ф(г1,г2)dгldг2
15 2 4^2-да-да(Г^г1 )(г2- ,) .
Тогда
F2 (щ,ш2 ) = ехр(1п F2 (ш-,ш2 )) = еФ еФ е_ Ф е_ Ф .
По формулам Сохоцкого
ф+—+ф—+ = 1 (—ф+^ 2ф) ,
с m-__L Г г ф(г1,г2)dhdh 12 ^2i-l(r1- t1 )(Г2- *2)
Следовательно, если —ф + ^ 2Ф = 0, то
F2 (^,^2) = еф еф
и решение краевой задачи (2.6) равно:
G++ (®j, ®2 ) = exp {ф++ (®j ,®2 )},
(2.7)
G (®j,®2) = ехр!Ф (®j,®2)}
После вычисления обратного преобразования Фурье получаем семейство решений исходной задачи. Единственное решение выбирается из условия выполнения равенства (0.2).
§3. Приближенные методы.
В параграфах 1-2 были изложены методы одновременного восстановления входного сигнала и аппаратной функции. Непосредственное применение предложенных формул возможно лишь в исключительных случаях, так как, как правило, преобразование Фурье и сингулярные интегралы в явном виде не вычисляются. В связи с этим возникает необходимость в разработке приближенных методов.
Одним из таких методов является вычисление интегралов (1.7)-(1.9) по квадратурным формулам вычисления преобразования Фурье [7], и вычисление сингулярных и бисингулярных интегралов в формулах
где
Ф++(Т1 ,Т2 ) = 1 [Ф (^Т2 )+ - | Ф^2 ^Ту +
4 7Т1 Т Тл
у —да 1 1
н± гФ^^т—_! г г______________фТг)_______аТ ат
^-да Т2 — Т2 ^2—I —1(Т1 —Т1) (Т21 — Т2 ) 1 2
(
^——! ^ 1 ^ N 1 да Ф(Т1,Т2),
Ф (Т1, Т2 )= . Ф(Т1, Т2) . I " ТТ1
4 у ™ —да т1 — т1
_! да Ф(Т1,Т2 )аТ — _! 7 да Ф(т1Т ) ТТ йт ™ —да Т1-Т2 ^2 —да—1(т1 — Т1 )(ти — Т2) 1 2,
из правых частей формул (2.7) по отимальным алгоритмам, предложенным в [8].
Более предпочтительным является итерационный метод, который обладает также свойством фильтрации. Вначале рассмотрим итерационный метод решения краевой задачи
1пО+ (ю) + 1п0~ (ю) = 1п^ (ю)
Воспользовавшись формулами Сохоцкого-Племеля, это уравнение представляется в виде
и(ю) = -TlU- (ю) + V(ю),
где
и +(м>) = 1п О +(м>), и— (^) = 1п(^),
V (^) = 1п (^).
Еще раз воспользовавшись формулами Сохоцкого-Племеля, получаем уравнение
(ю)
и(ю) =1 и(ю)—— Г и(т)Тт + V(ю),
'' ' ^ 7 2 ж? 1 т — ю v 7
—да
которые будем решать методом простой итера ип+1 (ю) =1 и„(ю)^т^ да 7^ + V(ю), (3.2)
... - (ю)
2т -ои
п = 0,1,....
Доказательство сходимости итерационного процесса 3.2 приводится в приложении 1.
Рассмотрим теперь итерационный метод решения уравнения (2.5).
Прологарифмировав уравнение (2.5), имеем
1п О++(ю1?ю2) + 1п О~~ (ю,ю)= 1п (юг,ю2) (3.3)
Введем функции
и++ (юг,ю2)=1пО++ (ю,ю), и~~ =,ю2)=1п = (ю,ю).
Представим уравнение (3.3) в виде
и++ (ю15ю2) + и~~ (юг,ю2)= 1п (юг,ю2)
Воспользовавшись формулами Сохоцкого [4] (в частности, используется формула и++— и+ — и += и ), предыдущее уравнение представим в виде
и (ю1?ю2) = 1п (ю,ю)—и +(ю,ю) — и_+(ю,ю).
Еще раз воспользовавшись формулами Сохоцкого [4], предыдущее уравнение можно записать в виде
1 ч 1
и(ю1?ю2) = 1п(юг,ю2) + 2 и(ю,ю)— ^(^и)(ю,ю),
где (^ )(»„»2 ) = ЧП и (Т"Т-)‘,Г,Г\ .
^ —1—1 (Т1 — ю1 )(Т2 — ю2 )
Функция и(ю1?ю2) вычисляется по итерационному процессу
ип+1 (ю1,ю2 ) = Кип (ю1,ю2 ) + (! — Кп )(ъ ^2 (ю1,ю2) + | (ю1, ю2 )- 1 (5цип )(ю!,ю2 (3-4)
п = 0,1,..., где 0<^<К <Р < 1.
Сходимость итерационного метода (3.4) к решению уравнения (3.3) доказана в приложении 2.
При практической реализации вычислительной схемы (3.4) проводится дискретизация и встречающиеся сингулярные интегралы вычисляются по оптимальным кубатурным формулам.
Рассмотрим модельный пример, иллюстрирующий эффективность предложенных алгоритмов. Дана модель, описываемая уравнением (2.1), где
а = 1,Т = (тх ,Т2 ),
И(т15Т2) = р(т,Г) при Т е[—1,1]2.
п2
И(Т1, Т2 ) = 0 при Т £[—1,1]2
х(^,Т2) = 0, при Т 0[—1,1]2, х (А, Тг ) = 1 при 0 < ^, Т2 < 1, х(^,Т2) = 2 при —1 < ^ < 0,0 < /2 < 1, х(Т\,Т2) = 1 при —1 <Ц < 0, — 1 < /2 < 0, х (А, Тг ) = —1 при —1 < ^ ,Т2 < 0,
Р(т, г) расстояние от точки Т до границы Г области [—1,1]2 .Функция
/(^,і2) = І І И(тх —іх,т2 —і2)• х(тх,т2)йтхйт2
—от —от
в результате ряда преобразований сводится к / (іх,12 ) = — (Х1 (^2,^2 ) + Х2 (^1,^2 ) + Х3 (^1,^2 ) + Х4 (^1, *2 )),
1 .у х 2
где ^1 (*1,ґ2) = |ії,у | х(у + *1,у + ґ2уіг , Х2 (^,*2) = || х(у + ^,г + /2)ф,
0 —у 0 —2
0 — у 0 — 2
X3 (^2) =| dу | х (у + г + У%2, X4 (ґ1,ґ2) =| <Лг | х (у + Ґ1, г + Ґ2 )/у.
—1 у —1 2
Требуется по функции / (іх ,і2) и условию
от от
| | |х(ї1,ґ2)|dt1dt2 = 5
—от —от
восстановить функции х (^ ,*2) и И (tх, *2 ).
Итерационный метод (3.4) позволил восстановить аппаратную функцию и входной сигнал с точностью 10—3.
Приложение 1. Докажем в метрике пространства Ь2 (—от, от) сходимость итерационного процесса (3.2). Известно (см.[ 9]), что
lU (t я
— J ШГ
mi J r — t
—ОТ
Следовательно
IK 1 (t) — U (t )||< 1|| Un (t) — Un_ x (t )|| + +1| U (t)—U„—1 (t )||< 3| |u (t)—U„—1 (t )| |.
Из этого неравенства и теоремы Банаха о сжатых отображения [10] следует сходимость итерационного
процесса (3.2) к решению U (t) уравнения
U * (t )=іU -(t)+v (t)
и справедливость оценки
\n
||и* (Т)— ип (Т)||< А(ЗМ)п.
Приложение 2. Докажем сходимость итерационного процесса (3.4) в метрике пространства Ь2. При этом
нам понадобятся следующие утверждения. Лемма. Справедливо равенство
11^12 ^)|| = |И.
Доказательство. Обозначим через Vtf преобразование Фурье функции /(тх,Т2) по переменной по переменной Т . Известно [9], что
V
1 Г^(г)
ui J Г — t
V —да
J------dr = F (rn)sgnrn
Кроме того, известно, что
\VУтf (ТТ) || = || f (т,т) 11.
Следовательно,
11^12 (и) || = \\К1К2 Я12 (и) || = || и(ю1,ю2>ёпю1 = II И(ю1,ю2) II = ||и(Т1,Т2) ||
т.е. Ц^иЦ =1.
Лемма доказана.
Теорема [11]. Пусть X — рефлексивное банахово пространство, A — линейный оператор с нормой И=1. Тогда итерационный процесс
xn+1 =KXn + (! — Лп) (Axn + f) ,
n = 0,1,..., где 0 <ы<Л<Р< 1 сходится к решению x уравнения x = Ax + f, если последнее имеет решение.
Нетрудно видеть, что оператор AU = — U S,U удовлетворяет условиям предыдущей теоремы, т.к. из
2 2 12
леммы следует, что \\AU\\=U. Отсюда следует сходимость итерационного процесса (3.4).
ЛИТЕРАТУРА
1. Сондхи М.М. // ТИИЭР - 1972.-N 7.-С.10 8-112.
2. Andrews Н.С. // Computer - 1974. N 5.-P. 36-39.
3. Бакут П.А., Свиридов К.Н., Сидельников В.Н., Устинов Н.Д. //Оптика и спектроскопия. 1982.-
Т.53. N 1.- С. 163-166.
4. Гахов Ф.Д. Краевые задачи.-М.: ГИФМЛ.-1963.-640с.
5. Гахов Ф.Д., Черский Ю.И. Уравнение типа свертки.-М.:Наука.-197 8.-2 9 6с.
6. Какичев В.А. Методы решения некоторых краевых задач для аналитических функций двух комплексных переменных. Тюмень.:Тюменский госуниверсистет. 1978.-124с.
7. Задирака В.К. Теория вычислений преобразования Фурье. Киев: Наукова думка. 1983-216с.
8. Бойков И.В.Оптимальные по точности алгоритмы приближенного вычисления сингулярных интегралов.-
Саратов: Изд-во Саратов. гос. ун-та 1983. 210с.
9. Иванов В.В. Теория приближенных методов и ее применение к численному решению сингулярных интегральных уравнений.-Киев. Наукова думка. 1968.-278
10. Люстерник Л.А., Соболев В.И. Элементы функционального анализа. М.: Наука, 1965. 540с.
11. Обломская Л.И. // ЖВМ и МФ 1968, т.8, N 2. c. 417-426.