Решение двухмерной коллимационной задачи рассеяния рентгеновских лучей с использованием нестандартных интегральных уравнений Текст научной статьи по специальности «Математика»

РЕШЕНИЕ ДВУХМЕРНОЙ КОЛЛИМАЦИОННОЙ ЗАДАЧИ РАССЕЯНИЯ РЕНТГЕНОВСКИХ ЛУЧЕЙ С ИСПОЛЬЗОВАНИЕМ НЕСТАНДАРТНЫХ ИНТЕГРАЛЬНЫХ УРАВНЕНИЙ Д.Д. Захаров, В.С. Сизиков, А.В. Смирнов, Б.А. Федоров

Сформулирована двухмерная задача о рентгеновском анизотропном рассеянии на объекте с использованием щелевого коллиматора и детектора. Выведено двухмерное сингулярное интегральное уравнение, связывающее истинную индикатрису рассеяния 3 с измеренной интенсивностью I. Уравнение приведено к нестандартной форме без сингулярности. Для его решения предложен итерационный метод. Дан численный пример.

Введение

Учет коллимационных искажений при малоугловом изотропном рассеянии рентгеновских лучей является обычной практикой во многих лабораториях [1-4]. При изотропном рассеянии на объекте относительная простота одномерного интегрального уравнения, связывающего экспериментальную I и искомую 3 (свободную от коллимационных эффектов) интенсивности рассеяния, позволила развить ряд методов его решения [1-4]. Эта простота объясняется тем, что интенсивость рассеяния не зависит от ориентации рассеивающего объекта относительно первичного пучка и, следовательно, является функцией лишь одного параметра, пропорционального при малоугловом рассеянии углу рассеяния.

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

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

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

Постановка задачи

Для повышения светосилы установки для изучения рассеяния рентгеновских лучей на образце часто используют щелевые коллиматоры и детекторы (рис. 1), что в отличие от точечных источников и приемников приводит к искажению («размытию») кривой рассеяния [1-5]. Чтобы восстановить истинное угловое распределение интенсивности рентгеновского рассеяния 3, которое было бы при точечном источнике и приемнике, нужно произвести коллимационный пересчет экспериментальной («размытой») интенсивности I [1, 3, 5-7].

Данная задача может быть одномерной, когда образец обладает изотропными свойствами. В этом случае обозначим через 3 (д) (или 3 (0) ) «истинную», неискаженную интенсивность (индикатрису) рассеяния рентгеновских лучей. Здесь д = (4тс/Х)б1п (0/2), где 0 - угол рассеяния, а X - длина волны рентгеновских лучей.

Задача может быть и двухмерной, когда свойства образца зависят от угла ф поворота оси симметрии образца относительно коллимационной щели (рис. 1). В этом случае искомая индикатриса рассеяния будет функцией двух переменных: 3 = 3 (д, ф) (или 3 = 3 (0, ф)).

Принципиальная схема хода лучей в рентгеновской камере представлена на рис. 1. Полагаем, что 0 << 1 (рассматривается малоугловое рентгеновское рассеяние). Параллельные рентгеновские лучи, формируемые коллимационной щелью длиной 2А1, падают на образец и рассеянные лучи, прошедшие через образец, регистрируются щелью детектора длиной 2А2, расположенной под углом 0. Использование щелей повышает светосилу установки, но регистрируемая щелью детектора интенсивность I = I(д, ф)

будет отличаться от 3 = 3(д, ф) , причем тем сильнее, чем больше длины щелей 2А1 и

2А2.

Ставится задача: восстановить математическим путем истинную индикатрису рассеяния 3 по измеренной функции рассеяния I. В одномерном случае такая задача уже рассматривалась в работах [1, 3, 5-7] и др., поэтому в данной работе она рассматривается кратко. Что же касается двухмерной задачи, то она рассматривается впервые в данной работе.

Одномерная задача

Пусть Б - область (прямоугольник) рассеяния на образце, Ь - расстояние от образца до плоскости регистрации, ё - ширина щели коллиматора, ё' - ширина щели детектора, g(х) - распределение интенсивности вдоль щели коллиматора, и(х') - распределение чувствительности вдоль щели детектора. Тогда суммарная интенсивность, регистрируемая всей щелью детектора при некотором угле 0 и при ё = ё', равна двойному интегралу по области Б и по щели детектора (ср. [5]):

< --(L) 1

2

J J(0) u(x') dx'

-Л2

g (x) dx,

(1)

где 0 - угол рассеяния от элемента dx области D в элемент dx' щели детектора.

Если распределения g (x) и u(x') равномерные, т.е. g (x) = const и u (x') = const.

то соотношение (1) может быть преобразовано в выражение (ср. [5], формула (21)):

/ (_ч

/ a < q < b, (2)

I (q) = 2 J J ((

2 + x2 )dx,

где a = qmm ^ 0, b = q

2 J J (

q2 + x2

> 0, или

dx = I (q), a < q < b.

(3)

Соотношение (3) является интегральным уравнением I рода типа Вольтерра, записанным в нестандартной форме, поскольку в нем нет ядра в явном виде, искомая функция J зависит не от одного внутреннего аргумента х, а от комбинации аргументов

■\]/2 + х2 и верхний предел интегрирования является некоторой функцией q. Теория и методы решения подобных уравнений практически не разработаны. Можно лишь отметить работы физического характера [1, 3] и работы [5-7] с математическим уклоном.

/2 2

Если ввести новую переменную ^ = ^/ + х , то уравнение (3) примет вид: ь

4

i

s2 - q2

J (s) ds = I(q), a < q < b.

(4)

Уравнение (4) является хорошо известным интегральным уравнением Абеля, принадлежащим к уравнениям типа Вольтерра I рода [8, с. 107], [9, с. 97]. Специфика уравнения (4) состоит в том, что оно является сингулярным: ядро зД/з2 - /2 обращается в бесконечность при ^ = / и это создает определенные трудности при его численном решении [10].

Уравнение (4) имеет аналитическое решение [5], [9, с. 98], [10] , _ 1'(/)

п

J (s) = - nj

П j

4q

2 2 2 - s2

dq,

a < s < b .

(5)

Однако вычисления по формуле (5) вызывают следующие осложнения. Во-первых, в (5) входит производная I'(/) от экспериментальной, а значит, зашумленной функции 1(/). Дифференцирование же зашумленной функции является некорректной (сильно неустойчивой) задачей и требует применения специальных устойчивых методов (регуляризации, фильтрации, сплайн-аппроксимации и т.д.) [8, 9, 11, 12]. Во-вторых, интеграл в (5) является сингулярным и для своего численного вычисления требует, например, использования обобщенных квадратурных формул [10]. От сингулярности в (5) можно избавиться, если выполнить интегрирование по частям. В этом случае [5], [9, с. 99]

1 ^(/) dq, а < з < Ь . (6)

J (s)=- JV^ dq

nJ dq

q dq

Однако решение (6) требует двукратного численного дифференцирования функции I(/). Чтобы преодолеть отмеченные трудности, в работе [10] предложен метод численного решения уравнения (4) на основе обобщенной квадратурной формулы, а

также метод численного вычисления интеграла в выражении (5) по обобщенной квадратурной формуле, учитывающей его сингулярность. А в работе [5] предлагается решать непосредственно нестандартное интегральное уравнение (3), причем решать его методом итераций.

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

Двухмерная задача

Пусть рассеивающие свойства образца зависят от ф - угла поворота оси симметрии образца, т.е. искомая индикатриса рассеяния является функцией двух переменных: 3 = 3(0, ф) или 3 = 3(д, ф). В этом случае суммарная интенсивность, регистрируемая щелью детектора, равна (ср. (1))

\2 4.1 ! г

g(х) ёх, (7)

' <0,')=(L П

L,

-Ai

J J(0, ф) u (x') dx'

где ф = ф + а (см. рис. 1). При u(x') = const и g (x) = const выражение (7) может быть приведено к виду (ср. (2)):

JbW _

1 (q ф) = J J (q2 + x2 , I arctg (x/q) + ф ) a ^ q ^ b, Фшт ^ ф ^ Фтах. (8)

Выражение (8) является интегральным уравнением относительно J при известной (измеренной) I, записанным в нестандартной форме. Однако выражение (8) является

неудобным для численной реализации, так как в нем аргумент sjq2 + x2 изменяется от b через q до b, т.е. немонотонно. Чтобы устранить этот недостаток, запишем (8) в виде: 1 (^ Ф) =

о __Vb2-q2 _

J J (q2 + x2 ,| arctg (x/ q) + ф|)) + J J (q2 + x2 ,| arctg (x/ q) + ф|)) =

-Vb^

л/Ь2-^^ __ТЪ2-? _

= | 3(д2 + х2 ,| arctg(х/д)-ф |)ёХ + | 3(д2 + х2,| аг^(х/д) + ф |)ёХ, о о

а < д < Ъ, фШ1П < ф < фтах . (9)

Метод итераций решения уравнения (9)

Для решения уравнения (9) предлагается, как и для решения уравнения (3) (см. [5]), использовать метод итераций Фридмана [8]:

лЪ-7

ji ф) = ji-i(q ф)+v

-q2

I(q, ф) - J Ji-1 (q2 + x2 , | arctg (x/q) - ф |)dx -

о

(10)

ib . _

J Ji-1 (q2 + x2 , | arctg (x/q) + ф |) dx

I = 1,2,3.....

о

где 0 < V < 211Щ. Норму оператора Ц можно оценить как (ср. [5])

||4 « 2^1 Ь2 - /2 < 2Ь , откуда

0 < V < 1/ Ь . (11)

Доказано [8, с. 272], что процесс итераций (10) сходится к точному решению в случае точной !(/,ф) при любом начальном приближении J0(q,ф) и при выполнении (11). Если же 1(/,ф) имеет погрешности, то процесс (10) расходится из-за некорректности задачи решения уравнения (9). В этом случае для выбора числа итераций можно воспользоваться известным правилом останова по невязке [8, с. 273-274], в которых число итераций согласуется с погрешностью измерений 5. Однако в данной работе мы касаться этого правила пока не будем.

Численный алгоритм

Численный алгоритм реализации схемы (10) является двухмерным обобщением итерационной схемы решения одномерного уравнения (2) или (3) согласно [5].

Введем равномерные совпадающие сетки узлов дискретизации по q и x: a, a + h1,K, b, где h1 = Лq = Л* = const - шаг дискретизации, или q( = a + h1 • i,

i = 0,1,..., M, xk = a + hb k, k = 0,1,..., M . Введем также равномерную сетку узлов по ф: Ф7- = ф0 + h2 • j, где h2 = Лф = const - шаг дискретизации по ф, j = 0,1,...,N .

Каждый из интегралов в (10) заменяем конечной суммой по формуле левых прямоугольников. При этом значения Ji-1 (л]q2 + x2, | arctg (x/q)+ ф |) , не попадающие в

узлы дискретизации, заменяем на Ji J nnd

Vq2 + x2

, nnd[ arctg(x/q)+ ф |] I, где через

nnd[-] обозначим значение ближайшего узла дискретизации (nearest node of discretization). В результате схема (10) в дискретном виде (при J0(q, ф) = 0) запишется как

J0(qi, ф j ) = 0

Ji(qi,фj) = Ji-1(qi,фj)+v|I(qi,фj)-h1 Z Ji-i^nnd Vq?+x2 ,nndlarctg(xkМ)-фj J

k=0

- h1 ZJi-1(nndl,nnd[arctg(x^qi )+фj |]U, i = 0,M, j = 0, N, i = 1,2,3, к

k=0

(12)

Разработана программа 1ТБЯ2 на языке ТигЬоС3 для решения уравнения (9) методом итераций согласно (12). С помощью программы 1ТЕЯ2 решен

Пример

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

(ф-п/ 2 )2

ехр -^^ , (13)

J ф) =

Г Л2 ' a \

\ q.

2а

3 3

где о = тс/5 ; параметры сетки узлов по д: а = дт1п = 5.5 • 10- , Ь = дтах = 16.5 -10" , шаг

Ы = Ад = 0.25 • 10-3, число шагов по д равно М = 44; параметры сетки узлов по ф: Фш;п = 0, фmax = П2, шаг Ь2 = Аф = п/90, число шагов по ф равно N = 45 ; 1/Ь = 60.6 (см. (11)); множитель V = 30 . На рис. 2 представлено точное решение 3 (д, ф) и измеренная функция I (д, ф) (найденная путем численного интегрирования формулы (9)).

Для характеристики точности итерационного процесса в зависимости от номера итерации I введем следующие функции: относительная погрешность 1-й итерации решения У/ (д, ф) по отношению к точному решению 3 (д, ф):

II Ф) -3Ф) IIь2

s = s отн l

II7

и относительная разность между l-й и (1-1)-й итерациями:

а = аотнl -•

ii 7i (q, ф) - 7i-i(q, ф)н l2

\Ji (q, ф)|Ь

Рис. 2. Точное решение J (q, ф) (а) и измеренная функция I (q, ф) (б)

На рис. 3 под цифрой 1 приведены полученные зависимости sотнi и аотнi. Имеем: s - min - 0.06730 при l - 29, а - min - 2.761 -10-3 при l - 25 .

На рис. 4 приведено решение J29(q, ф), соответствующее минимуму sотнl.

Далее с помощью датчика случайных чисел RNDAN [9, с. 153] к значениям I (q, ф) были добавлены нормальные погрешности с нулевым математическим ожиданием и среднеквадратическим отклонением аi - 0.0001, что соответствует относительной погрешности || AI ||/1 - 0.03005 « 3%. На рис. 5а приведена зашумленная I (q, ф). На рис. 3 под цифрой 2 приведены зависимости sотнl и аотнl. Имеем:

s - min - 0.1660 при l - 9, а - min - 0.01642 при l - 15 . На рис. 5б приведено решение J9(q, ф), соответствующее минимуму s^^l.

2

2

Рис. 3. Функции s = вотнi (непрерывные кривые) и а = аотнi (пунктирные кривые)

э б

Рис. 5. Зашумленная I (q,ф) (а) и решение Jg(q, ф) (б)

Для улучшения решений Ji (q, ф) выполнено сглаживание функции I (q, ф) окном 3 х 3 . На рис. 6а приведена сглаженная I (q, ф), обозначенная как S(q, ф) . На рис. 3 под цифрой 3 приведены зависимости Sothi и аотнi. Имеем: s = min = 0.1125 при i = 13 ,

а- min - 9.331-10 при l - 16. На рис. 6б приведено решение вующее минимуму sотнl.

J13(q, ф), соответст-

Рис. 6. Сглаженная I (q, ф) = S(q, ф) (а) и решение J13(q, ф) (б)

Для дальнейшего улучшения решений Jl (q, ф) выполнено сглаживание не только функции I (q, ф) окном 3 х 3 , но и сглаживание решения Jl (q, ф) в каждой итерации (также окном 3 х 3 ). На рис. 3 под цифрой 4 приведены зависимости sотнl и аотнl.

Имеем: s - min - 0.09737 при l - 28, а - min « 8-10_4 при l > 100 . На рис. 7 приведено решение J28(q, ф), соответствующее минимуму sотнl.

Рис. 7. Решение ^28(д, ф) при сглаживании I (д, ф) и у (д, ф)

Отметим, что результаты, приведенные на рис. 4—7, получены при У)(д, ф) = 0, т.е. при нулевом начальном приближении. Можно еще улучшить решения 31 (д, ф) за счет изменения начального приближения. Было использовано следующее начальное приближение:

J (q, ф) =

f 6 -10-3 ^ Г (ф-П2 -П90)2

- exp - —--2—

t q J L 2(п/4)2

весьма близкое к искомому решению (13). На рис. 3 под цифрой 5 приведены зависи-

—3

мости вотнi и аотнi. Имеем: s = min = 0.08013 при l = 3, а = min = 2.266 -10 при l = 17 . На рис. 8 приведено решение J3(q, ф), соответствующее минимуму 8отнi.

(14)

Рис. 8. Решение /з(д, ф) при сглаживании I (д, ф) и / (д, ф) и начальном приближении (14)

Заключение

Сформулирована двухмерная задача о анизотропном рентгеновском рассеянии на объекте, когда измеренная интенсивность I и искомая индикатриса рассеяния / зависят от угла рассеяния 9 и от угла поворота оси симметрии образца ф. Для восстановления / по экспериментальной I выведено двухмерное сингулярное интегральное уравнение. Уравнение приведено к нестандартной форме без сингулярностей. Для его решения предложен итерационный метод. Приведен численный пример.

Результаты решения численного примера, а также ряда других (неприведенных в работе) примеров показывают следующее. Применение итерационной схемы Фридмана (10)-(12) для решения неклассического двухмерного интегрального уравнения (9) является весьма эффективным. Выбор оптимального числа итераций п8 0рг, соответствующего минимуму функции в = 8отн/, на практике невозможен, так как эта функция

включает в себя точное (неизвестное) решение / (д, ф). На практике можно использовать лишь функцию а = аотн/. Как показывает рис. 3 и решение других примеров, число итераций па ^ в среднем несколько больше п8 0рг, тем не менее, па ^ можно использовать для оценки п8 ^. Точность итерационного решения / (д, ф) повышается

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

Литература

1. Guinier A., Fournet G. Small-angle scattering of X-rays. - NY: Wiley, 1955.

2. Федоров Б.А. Учет коллимационных искажений при малоугловом рассеянии рентгеновых лучей. Поправка на высоту щелей // Кристаллография. 1968. Т. 13, № 5. С. 763-769.

3. Schelten J., Hossfeld F. Application of spline functions to the correction of resolution errors in small-angle scattering // J. Appl. Cryst. 1971. Vol. 4. P. 210-223.

4. Мельничук А.П., Прищепенок О.Б., Смирнов А.В., Федоров Б.А. Прецизионная юстировка камеры Краткого и программа первичной обработки данных рентгеновского малоуглового рассеяния // Изв. вузов. Приборостроение. 2002. Т. 45, № 7. С. 4854.

5. Сизиков В.С., Смирнов А.В., Федоров Б.А. Решение одномерной коллимационной задачи оценки рентгеновского изотропного рассеяния излучения методом итераций // Изв. вузов. Приборостроение. 2005. Т. 48, № 10. С. 44-52.

6. Добровольский В.А., Гориловский Д.А., Сизиков В.С., Смирнов А.В., Федоров Б.А. Модификация метода квадратур численного решения интегрального уравнения Абеля с учетом его сингулярности // В сб. научн. статей «Современные технологии» под ред. С. А. Козлова. СПб: СПбГИТМО(ТУ), 2001. С. 121-126.

7. Гориловский Д.А., Добровольский В.А., Сизиков В.С., Смирнов А.В., Федоров Б.А. Решение интегрального уравнения Абеля в нестандартной форме методом итераций // В сб. научн. статей «Современные технологии» под ред. С. А. Козлова. СПб: СПбГИТМО(ТУ), 2001. С. 127-134.

8. Верлань А.Ф., Сизиков В.С. Интегральные уравнения: методы, алгоритмы, программы. Киев: Наук. думка, 1986.

9. Сизиков В.С. Математические методы обработки результатов измерений. СПб: Политехника, 2001.

10. Сизиков В.С., Смирнов А.В., Федоров Б.А. Численное решение сингулярного интегрального уравнения Абеля обобщенным методом квадратур // Изв. вузов. Математика. 2004. № 8(507). С. 62-70.

11. Петров Ю.П., Сизиков В.С. Корректные, некорректные и промежуточные задачи с приложениями. СПб: Политехника, 2003.

12. Petrov Yu.P., Sizikov V.S. Well-Posed, Ill-Posed, and Intermediate Problems with Applications. Leiden-Boston: Brill Acad. Publishers, VSP et al., 2005.