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

Бойков И.В., Мойко Н.В. ОБ ОДНОМ МЕТОДЕ ВОССТАНОВЛЕНИЯ ИМПУЛЬСНОЙ ФУНКЦИИ И ИСКАЖЕННОГО АТМОСФЕРОЙ КОРОТКОЭКСПОЗИЦИОННОГО ИЗОБРАЖЕНИЯ

Введение.

Изображение объекта, наблюдаемого через турбулентную атмосферу, является сильно искаженным. Для восстановления изображения существуют различные методы фильтрации [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.