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

ПРИМЕНЕНИЕ МЕТОДА ИТЕРАЦИЙ ДЛЯ ПОВЫШЕНИЯ ТОЧНОСТИ ОПРЕДЕЛЕНИЯ ИНДИКАТРИСЫ РАССЕЯНИЯ РЕНТГЕНОВСКИХ ЛУЧЕЙ Д.Д. Захаров

Научные руководители - доктор технических наук, профессор В.С. Сизиков, кандидат физико-математических наук, доцент А.В. Смирнов, доктор физико-математических наук, профессор Б.А. Федоров

Сформулирована двухмерная задача о рентгеновском рассеянии на анизотропном объекте при использовании щелевого коллиматора и детектора. Для восстановления истинного углового распределения интенсивности рентгеновского рассеяния 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) ) «истинную», неискаженную ин-

тенсивность (индикатрису) рассеяния рентгеновских лучей. Здесь q = (4тс/^)вт(0/2), где 0 - угол рассеяния, а X - длина волны рентгеновских лучей.

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

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

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

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

Рис. 1. Схема хода лучей в рентгеновской камере

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

сивность, регистрируемая всей щелью детектора при некотором угле 0 и при d = d', равна двойному интегралу по области Б и по щели детектора (ср. [5]):

А,

Г A

=(£ I

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

g( x) dx, (1)

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

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

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

_

I(q) = 2 J J (q2 + x2 )dx, a < q < b, (2)

0

где a = qmin ^ 0, b = qmax > 0, или

Vb^-q2 _

2 I J (q2 + x2 )dx = I(q), a < q < b. (3)

0

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

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

V2 2

q + x , то уравнение (3) примет вид:

b

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

Js2 - q2

4

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

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

3(я) = -1 [ / (д) dq, а < я < Ь . (5)

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

b ' 1 dI (q)

j (s)=- Ivq^ d

nJ dq

q dq

dq, a < s < b . (6)

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

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

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

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

Пусть рассеивающие свойства образца зависят от ф - угла поворота оси симметрии образца, т.е. искомая индикатриса рассеяния является функцией двух переменных: J = J(0, ф) или J = J(д, ф) . В этом случае суммарная интенсивность, регистрируемая щелью детектора, равна (ср. (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 (^ Ф) = J J ((2 + x2, I arctg (x/q) + ф ) a < q < Ъ, pmin < ф < фтах. (8) -jtw

Выражение (8) является интегральным уравнением I рода относительно J при известной (измеренной) I, записанным в нестандартной форме. При этом выражение (8) неудобно для численной реализации, так как в нем аргумент л]q2 + x2 изменяется от Ъ через q до Ъ, т.е. немонотонно. Чтобы устранить этот недостаток, запишем (8) в виде: 1 (^ Ф) =

о __Vb2-q2 _

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

~7 0

ТЪч

W __. . _

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

Vb^-q2 __V^v

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

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

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

jiф) = ji-ife ф)+v

-q2

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

о

(10)

4b _ _

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

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

где 0 < v < 2 /1|4||. Норму оператора можно оценить как (ср. [5]) 114 « 2^1 b2 - q2 < 2b ,

откуда

0 < v < 1/ b . (11)

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

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

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

Введем равномерные совпадающие сетки узлов дискретизации по q и х: a, a + h1,..., b, где h1 = Aq = Ax = const - шаг дискретизации, или q( = a + h1-i,

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

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

дискретизации, заменим на Ji_1 (nnd ^q2 + х2 ,nnd[ | arctg (x/q) + ф |]) , где через nnd[-]

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

J0( qi, ф j ) = 0,

Г M

Ji(qi,ф]) = Ji_1(qi,ф])+vjI(q,ф;)_h1^Ji_1 (nnd Vqf + xk ,nnd[|arctg(xk/q)_ф;|])

к=0

M

-h1^ Ji_1 (nnd V?

+ xk

nnd

[| arctg(xk/qt) + ф] |])|, i = 0,M, j = 0,N, i = 1,2,3,.

к=0

(12)

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

Пример

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

)2 ' (13)

J ф) =

Г Л2 ' a \

\ q,

exp

(ф_П 2 )2

2а

где а = п/5; параметры сетки узлов по q: a = qmin = 5.5 -10 , b = qmax = 16.5 -10 , шаг

h1 = Aq = 0.25 • 10 3, число шагов по q равно M = 44; параметры сетки узлов по ф:

фт1п = 0, фтах = V2, шаг Ъ2 = Дф = п/90, число шагов по ф равно N = 45 ; 1/Ъ = 60.6 (см. (11)); множитель V = 30. На рис. 2 представлено точное решение J(д,ф) .

и-1

Рис. 2. Точное решение J (д, ф)

На рис. 3 приведена измеренная функция I (д, ф).

Рис. 3. Измеренная функция I (д, ф)

На рис. 4 приведены следующие функции: относительная погрешность 1-й итерации решения J^ (д, ф) по отношению к точному решению J (д, ф) (непрерывная кривая)

II Jl(д, ф) - Уф)

8 = 8

отн I

и (д, ф)1

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

ь

2

а = а

отн l

II Ji (q ф) - Jiф) llz2 II Ji (q, фЖ,

0.9

0.3

0.2

1

|\

1 \ \ N

\ ч ,

0 10 20 30 40 50 60 70 80 90 l 100

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

кривая)

В данном примере получен следующий результат: s = min = 0.0673 при l = 29.

а = min = 2.761 • 10~3 при 1 = 25.

Рис. 5. Решение J29(q, ф)

На рис. 5 приведено решение ф) • Сравнение рис. 2 и 3 говорит о том, что численное решение ф) дает удовлетворительное приближение к точному решению 1 ^ ф).

0.8

0.7

0.6

0.5

0.4

0.1

0

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

Заключение

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

Литература

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

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

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. С. 48-54.

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

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

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

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

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

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