Введение.

В задачах радиографии информацию о внутренней структуре объектов получают путем просвечивания этих объектов проникающим излучением и анализа прошедшего излучения, которое состоит из первичной (нерассеянной) и рассеянной компонент. Традиционные методы и схемы радиометрической дефектоскопии (см., например, [1]) учитывают только первичную компоненту при анализе и интерпретации прошедшего излучения, при этом используются различные методики уменьшения влияния рассеянного излучения (нежесткая коллимация, дискриминация по фотопику и др.). Как правило, применение этих методик требует существенного увеличения мощности проникающего излучения и усложнения способов регистрации прошедшего сквозь объект излучения. Между тем, детальное знание распределения рассеянного излучения позволяет корректно учесть его при использовании радиографических методов, ориентированных на обработку первичного излучения, и дает возможность снизить экспозиционную дозу облучения при проведении экспериментальных работ. Кроме того, в некоторых случаях (например, при исследовании структуры труднодоступных образцов) рассеянное излучение является единственным источником информации об объекте при их радиографическом анализе.

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

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

Например, широко распространенный программный комплекс MCNP [2] (The N-particle Monte-Carlo code) дает возможность исследовать процессы формирования образов с учетом как первичного, так и рассеянного излучения и провести анализ влияния различных механизмов взаимодействия фотонов с веществом (фотопоглощения, когерентного и Комптоновского рассеяния, образования электрон-позитронных пар и т.д.). Однако, описание геометрии сложных многокомпонентных объектов (таких, например, как литье) является чрезвычайно сложным в рамках указанного программного продукта. Кроме того, специфика использованного в MCNP подхода к распараллеливанию накладывает обременительные ограничения на число фотонных историй при проведении вычислительных экспериментов на большинстве используемых 32-х разрядных процессорных системах, что может привести к слишком большой погрешности получаемых результатов.

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

Физическая модель.

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

В настоящей работе рассматривается диапазон энергий фотонов проникающего излучения до 1 МэВ (Рентгеновское излучение). Основными процессами взаимодействия квантов с веществом в этом диапазоне являются когерентное (Рэлеевское) и некогерентное (Комптоновское) рассеяние и фотоэлектрическое поглощение фотонов (рис.1). Эффекты связанности электронов в атоме учитываются путем введения форм-фактора для когерентного и функции рассеяния для некогерентного рассеяний [5,6]. Вторичные эффекты, такие как флуоресценция и перенос электронов, несущественны в рассматриваемом энергетическом диапазоне с точки зрения влияния на исследуемое распределение прошедшего сквозь объект излучения (см., например [2]).

Как известно [7,8], если рассматривать рассеяние фотона на электроне, свободном от связи с другими частицами, то дифференциальное сечение рассеяния неполяризованного фотона  зависит только от энергии и направления рассеиваемого фотона и выражается формулой Клейна-Нишины-Тамма. Предположение о свободном и покоящемся электроне справедливо только при условии, что импульс, передаваемый электрону, намного превышает импульс его первоначального движения в атоме , где  и  ‑ первоначальные длины волн соответственно фотона и электрона,  - угол рассеяния фотона. Наличие связей и первоначального движения электронов приводит к уменьшению вероятности комптоновского рассеяния. В этом случае сечение комптоновского рассеяния определяется комбинацией  и нерелятивистской функции некогерентного рассеяния Хартри-Фока  [5]:

(1)

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

Рассеяние фотонов на электронах становится когерентным при небольшой величине передаваемого импульса (Рэлеевское рассеяние). Такое рассеяние происходит без потери энергии квантом проникающего излучения. Дифференциальное сечение Рэлеевского рассеяния в случае учета связанности электронов является произведением Томпсоновского классического сечения  и квадрата релятивистского атомного форм-фактора Хартри-Фока  [6]:

,

(2)

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

На рис.3, 4 изображены графики функции,  и  для алюминия и железа ().

Подробно используемая физическая модель описана в [9].

Методика моделирования процессов рассеяния излучения.

Для исследования взаимодействия ионизирующего излучения с гетерогенными кусочно-однородными объектами и анализа закономерностей формирования рентгеновских изображений их внутренней структуры разработан расчетный аппарат на основе [3]. Этот аппарат включает в себя:

-       Поверхностно ориентированное описание многокомпонентных объектов сложной формы, основанное на точном определении границ гомогенных составляющих внутренней структуры объекта и задании соответствующих разделяющих замкнутых оболочек;

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

1. Гетерогенные среды с кусочно-однородной структурой могут быть адекватно описаны посредством точного задания разделяющих гомогенные части таких сред поверхностей. Такое поверхностно ориентированное описание кусочно-однородных объектов позволяет определять соответствующие структуры с использованием замкнутых оболочек. Для дискретизации замкнутых оболочек, описывающих гомогенные компоненты объекта, в описываемом алгоритме используется триангуляционная модель.

Необходимо отметить, что подход, основанный на описании граничных поверхностей, имеет важную особенность, которая должна быть учтена при разработке эффективного вычислительного алгоритма, а именно:

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

Отметим еще одну важную черту, характерную для рассматриваемого подхода. Количество рассматриваемых лучей может достигать десятков и сотен миллионов для каждой детекторной плоскости, а разделяющая оболочка после триангуляции может состоять из десятков тысяч треугольников для достаточно сложных структур. Такие величины параметров дискретной модели требуют разработки эффективного способа просмотра всех возможных пар «луч-треугольник» при поиске точек пересечения луча с поверхностью разделяющей оболочки.

Предложен следующий двухэтапный алгоритм, учитывающий указанные выше особенности рассматриваемого поверхностно ориентированного подхода:

-       сначала определяются все треугольники, которые пересекает данный луч. Этот «качественный» этап проводится с использованием рамки (ребер и вершин) треугольника, а факт пересечения лучом треугольника устанавливается с помощью специальной системы неравенств (см. [3]);

-       затем с использованием фасеток вычисляются координаты точек пересечения («количественный» этап) данного луча с найденными на первом этапе треугольниками.

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

2. Рассмотрим теперь основные особенности реализации метода Монте-Карло с использованием «оболочечного» описания объектов.

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

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

Сначала определяется оптическая толщина объекта в направлении движения фотона от начала звена до выхода из объекта, при этом используется разработанный в [3] алгоритм нахождения точек пересечения соответствующего луча с разделяющими оболочками.

После определения конца текущего звена траектории в случае, если фотон остался в объекте, проводится розыгрыш типа взаимодействия этого фотона с веществом той части объекта, в которой заканчивается рассматриваемое звено траектории.

Следующим этапом алгоритма является определение направление движения фотона после рассеяния (в случае, если в результате взаимодействия этот фотон не поглотился). Это направление вычисляется с помощью розыгрыша случайной величины – угла рассеяния кванта с использованием известных угловых распределений, соответствующих рассматриваемым типам рассеяния.

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

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

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

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

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

Отметим в заключение, что описанный в работе алгоритм реализован в виде параллельного программного комплекса с использованием интерфейса MPI, ориентированного на многопроцессорные вычислительные системы. Расчеты проводились на МВС-5000М.

Результаты моделирования.

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

Ниже обсуждаются примеры расчета трансформации рентгеновского излучения в объектах, дается феноменологическая интерпретация некоторых эффектов.

Рассмотрим вначале результаты сравнительных расчетов, проведенных с использованием описанного выше алгоритма и с помощью известного пакета MCNP. Для сравнения в расчетах использовался монохроматический источник рентгеновского излучения с энергией фотонов 100 кэВ, железная пластина 5х5х1 см. Результаты расчетов как прямого, так и рассеянного излучений (рис.5, 6) совпадают в рамках статистической погрешности расчетов.

Рис.5	Рис.6

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

Рис. 8                                                                              Рис. 9

Рассмотрим теперь влияние рассеянного излучения на качество рентгенографического изображения «дефекта» на примере, в котором железная пластина с полостью («дефектом») просвечивается излучением с энергией 80 кэВ (рис.7). В расчетах использовалась пластина 5х5х1 см с полостью посередине 0.8х0.8х0.4 см. Источник излучения расположен сверху пластины на расстоянии 100 см. Детекторная плоскость расположена вплотную к нижней грани объекта и имеет 300х300 детекторов размером 0.5х0.5 мм. На рис.8,9 приведены цифровые изображения рассматриваемого объекта в прямом (рис.8) и рассеянном (рис.9) излучении для случая, когда детекторы регистрируют интенсивность прошедшего сквозь пластину излучения. Эти иллюстрации показывают, что для рассматриваемой ситуации изображение полости (дефекта), формируемое прямым (нерассеянным) излучением практически отсутствует, в то время как в рассеянном излучении эта полость достаточно хорошо видна (светлые участки соответствуют большей интенсивности, темные – меньшей). Указанные отличия в изображении полости в прямом и рассеянном излучении объясняется следующим образом.

На рисунках 10 и 11 изображены графики зависимости интенсивности регистрируемого излучения  при  ( - координата детектора по горизонтали,  - по вертикали).

Рис. 10

Рис. 11

Из этих графиков видно, что относительное изменение интенсивности прямого излучения (в окрестности ), обусловленное наличием полости, мало заметно на фоне максимальной интенсивности (рис.10), в то время как аналогичное изменение в рассеянном излучении (рис.11) весьма велико.

Отметим также хорошо заметное подчеркивание границ объекта на изображении, формируемом рассеянным излучением (рис.9). Рассмотрим этот эффект более подробно.

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

В области II () объем рассеивающего вещества, дающий вклад в соответствующие детекторы, уменьшается с приближением к границе пластины до половины полусферы для точек, лежащих на границе. Наконец, в области III () вклад в детекторы вносит излучение, рассеянное в прямоугольной области  по всей толщине пластины, поскольку вылетевшие за пределы объекта фотоны не испытывают поглощения. Поэтому происходит увеличение интенсивности регистрируемого излучения. Затем эта величина достигает максимума и уменьшается в соответствии с угловым распределением рассеянных фотонов. Описанная выше закономерность формирования радиографического образа объекта в рассеянном излучении хорошо видна на рис.11 в окрестности края пластины (). Величина максимума и его положение зависит от соотношения величин  и . В случае  максимума интенсивности нет. Детально рассмотренные эффекты обсуждаются в [10].

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

Литература.

1.   Радиоизотопная дефектоскопия. М., Атомиздат, 1976.

2.   J.F.Briesmeister (ed.). MCNP - A General Monte Carlo N-Particle Transport Code. LANL Report LA-13709-M, Los Alamos, 2000.

3.   В.П.Загонов, М.Е.Жуковский, С.В.Подоляко, М.В.Скачков, Г.-Р.Тиллак, К.Беллон. Применение поверхностно ориентированного описания объектов для моделирования трансформации рентгеновского излучения в задачах вычислительной диагностики. М., Математическое моделирование. 2004, т.16, №5,.с.103-116.

4.   Proceedings of the 2^nd International Conference on Computer Methods and Inverse Problems in Nondestructive Testing and Diagnostics. Minsk, Belarus, DGZfP Berichtsband 64, 1998.

5.   J.H.Hubbell, W.J.Veigele, E.A.Briggs, R.T.Brown, D.T.Cromer, and R.J.Howerton. Atomic Form Factors, Incoherent Scattering Functions, and Photon Scattering Cross Sections. J. Phys. Chem. Ref. Data 4, 471-538 (1975); erratum in 6, 615-616 (1977).

6.   J.H.Hubbell, and Overbo. Relativistic Atomic Form Factors and Photon Coherent Scattering Cross Sections. J. Phys. Chem. Ref. Data 8, 69-105, (1979).

7.   У.Фано, Л.Спенсер, М.Бергер. Перенос гамма-излучения. Госатомиздат, Москва, 1963.

8.   Экспериментальная ядерная физика, под ред. Э. Сегре, Издательство «Иностранной литературы», Москва, 1955.

9.   С.В.Подоляко, Е.Г.Лукьянова. Численное моделирование трансформации рентгеновского излучения в объектах с учетом  влияния форм‑факторов  на угловое распределение фотонов. Препринт ИПМ им. М.В.Келдыша РАН №6 за 2004г.

10.             V.M.Artemiev, J.N.Gray, G.-R.Tillack. Scattering effect in radiation techniques – theory and experiment. Computer Methods and Inverse Problems in Nondestructive Testing and Diagnostics. The 2^nd International Conference Minsk, Belarus, DGZIP, Berichtsband 64, 323-330, 1998.