УДК 620.179.152
ПОВЫШЕНИЕ ПРОНИКАЮЩЕЙ СПОСОБНОСТИ СИСТЕМ ЦИФРОВОЙ
РАДИОГРАФИИ НА ОСНОВЕ АНАЛИЗА СИГНАЛОВ НИЗКОЙ
ИНТЕНСИВНОСТИ
© 2022 г. В.Ю. Жвырбля1,* С.П. Осипов1,**, Д.А. Седнев1,***
1Томский политехнический университет, Россия 634028 Томск, пр-т Ленина, 30
E-mail: * vzhrb3001@gmail.com; ** osip1809@rambler.ru; ***sednev@tpu.ru
Поступила в редакцию 16.05.2022; после доработки 30.05.2022
Принята к публикации 03.06.2022
Разработана математическая модель формирования высокоэнергетических цифровых радиографических изображе-
ний в условиях низких уровней цифровых сигналов. Модель предназначена для повышения эффективности анализиру-
емых систем применительно к контролю крупногабаритных объектов за счет цифрового суммирования изображений. На
базе модели разработан алгоритм и программа на MathCad, предназначенные для симуляции цифровых изображений
испытуемых объектов. На основе анализа результатов вычислительных и натурных экспериментов проиллюстрирована
принципиальная возможность повышения эффективности систем высокоэнергетической цифровой радиографии на
основе учета редких событий.
Ключевые слова: тормозное излучение, цифровая радиография, проникающая способность, разрядность АЦП, макси-
мальная энергия тормозного излучения, редкие события, имитационное моделирование.
DOI: 10.31857/S0130308222070041, EDN: BOEOUQ
ВВЕДЕНИЕ
Расширяющиеся потребности промышленности и науки в высококачественном контроле и
исследовании крупногабаритных объектов обусловливают необходимость совершенствования
методов и средств неразрушающих испытаний [1—5]. Применение систем цифровой радиогра-
фии (ЦР) и рентгеновской компьютерной томографии (КТ) ограничивается двумя основными
факторами. Первый фактор имеет исключительно физическую природу и обусловлен малой про-
зрачностью исследуемого объекта для тормозного излучения [6, 7]. Отмеченный вид прозрач-
ности объекта контроля (ОК) называется радиационной (радиоскопической) прозрачностью [8,
9]. Второй фактор связан с малым уровнем сигналов на выходе детекторов [10, 11], формирую-
щих изображения в ЦР или проекций в КТ. Существует ряд причин, приводящих к малому уров-
ню сигналов. Эти причины связаны с параметрами ОК, схемы сканирования, характеристиками
источника и детектора излучения. Основными параметрами источника излучения применитель-
но к обсуждаемому вопросу являются: для счетного режима регистрации фотонов — поток
частиц; для спектрометрического режима — распределение потока фотонов по энергии; для
интегрального режима — интенсивность излучения. Объект контроля характеризуется толщи-
ной, плотностью и эффективным атомным номером, впрочем, эти три параметра сводятся к двум
— массовой толщине и эффективному атомному номеру материала или к одному — толщине ОК
в длинах свободного пробега (д.с.п.) для заданной максимальной энергии тормозного излучения.
Под массовой толщиной понимается произведение толщины на плотность. Толщина ОК в д.с.п.
является наиболее удобным параметром для оценки предельных возможностей систем ЦР и КТ
применительно к контролю крупногабаритных объектов. В счетном и спектрометрическом
режиме пересчитываются фотоны, провзаимодействовавшие с чувствительным объемом детек-
тора, только в первом случае — гипотетически для всего диапазона изменения энергии зареги-
стрированного фотона, а во втором — фотоны распределяются по энергетическим окнам.
Интегральный режим регистрации характеризуется разрядностью аналого-цифрового преобра-
зователя (АЦП) и реальным максимальным значением цифрового сигнала (ЦС). Для всех режи-
мов регистрации малые значения сигналов связаны, прежде всего, с малым числом фотонов,
провзаимодействовавших с чувствительным элементом детектора. В связи с изложенным выше,
отмеченную ситуацию называют «фотонным голоданием» («photon starvation») [12, 13]. Для
интегрального режима регистрации помимо «фотонного голодания» к малым значениям ЦС при-
водят недостаточная разрядность АЦП и неполное использование диапазона изменения ЦС [11,
14—16]. Причем степень отмеченной недостаточности разрядности существенно зависит от
вида решаемой задачи испытаний, исследуемого объекта, заданного качества контроля (произ-
водительность, точность измерений, разрешение по эффективному атомному номеру в досмо-
40
В.Ю. Жвырбля, С.П. Осипов, Д.А. Седнев
тровом контроле и т.п.). В работах [11, 14, 17] подчеркивается необходимость согласования
диапазонов изменения аналоговых сигналов (АС) и ЦС и отмечается, что лучшего согласования
можно достичь в случае использования в детекторах, работающих в интегральном режиме реги-
страции, логарифмических усилителей, предшествующих АЦП. Причем для логарифмических
усилителей требуются АЦП с существенно меньшей разрядностью. Отметим, что «фотонное
голодание» устраняется или существенно снижается, если нет сопутствующих технических
ограничений, тремя основными способами: увеличением максимальной энергии тормозного
излучения; повышением мощности источника фотонов; увеличением времени измерения. Для
интегрального режима регистрации фотонов изменение параметров систем ЦР и КТ принесет
заметный эффект исключительно в случае изначального существенного несоответствия диапа-
зонов изменения АС и ЦС. Для интегрального режима регистрации необходимо согласовывать
диапазоны изменения АС и ЦС особенно в случае возможности существенного уменьшения
эффекта «фотонного голодания».
Если замена источника излучения на более интенсивный невозможна, а повышение энергии
тормозного излучения ограничено фиксированным уровнем, то актуален вопрос о применении
других подходов к повышению качества формируемых цифровых радиографических изображе-
ний или проекций. Для интегрального режима регистрации искомый подход может быть связан
с представлением оцифрованного АС случайной величиной и анализом ее для описанного выше
случая. Для проверки обоснованности отмеченного подхода необходима разработка математиче-
ской модели формирования радиографических изображений (проекций) применительно к испы-
танию крупногабаритных ОК с учетом редких событий. Редкие события возникают в случае
отличия ЦС от нуля при условии, что средние значения соответствующих АС могут быть близки
к аналоговому эквиваленту цифровой единицы, но меньше ее. Если не учитывать этот фактор,
то информативность полученных цифровых радиографических изображений будет нулевой.
Представляют интерес и условия наступления редких событий, которые обуславливают практи-
ческую применимость математической модели.
1. МАТЕМАТИЧЕСКАЯ МОДЕЛЬ ФОРМИРОВАНИЯ СИГНАЛОВ
При разработке математической модели формирования ЦС для систем ЦР с импульсными
источниками тормозного излучения возьмем за основу [18] с учетом особенностей, отмеченных во
введении.
1.1. Параметры математической модели
Система высокоэнергетической цифровой радиографии характеризуется параметрами, связан-
ными с ОК, источником и детектором тормозного излучения.
К параметрам математической модели относятся:
максимальная энергия тормозного излучения Emax, МэВ;
числовой энергетический спектр источника тормозного излучения f(E, Emax);
частота следования импульсов ν, Гц;
среднее количество фотонов, испускаемых источником и попадающих на фронтальную
поверхность чувствительного объема сцинтиллятора, N1;
количество импульсов для формирования цифрового сигнала n;
химическая формула вещества чувствительного элемента детектора Chd, которая в обобщен-
ном виде представляется выражением:
n
d
Chd
(
A
,
p
)
,
(1)
=
di
di
i=1
где nd — количество химических элементов, образующих вещество сцинтиллятора; Ad i, pd i — обо-
значение химического элемента, его количество атомов в молекуле;
толщина чувствительного объема детектора hd, см;
плотность материала чувствительного объема детектора ρd, г/см3;
разрядность АЦП kADC, бит;
степень заполнения диапазона изменения цифрового сигнала CADC;
максимальная толщина ОК в длинах свободного пробега P(Emax), д.с.п.;
минимальный уровень цифрового сигнала M.
Дефектоскопия
№ 7
2022
Повышение проникающей способности систем цифровой радиографии на основе анализа...
41
Параметрами материала чувствительного объема детектора являются векторы молярных масс
Md = (Md 1, Md 2,…, Md n) и атомных номеров элементов Zd = (Zd 1, Zd 2,…, Zd n), которые однозначно
соответствуют вектору Ad = (Ad 1, Ad 2,…, Ad n) [19].
1.2. Основные формулы
Выражение, связывающее параметры математической модели и величину цифрового сигнала
D, имеет вид:
n
N
1i
(P)
D(P)
=
C
E (P)+E
∆
,
(2)
∑
E
∑
j
bi
1
i
=1
j
=1
где [arg] — целая часть arg; CE — коэффициент преобразования энергии фотонов, поглощенной в
чувствительном объеме, в электрическую энергию; Ni(P), i = 1…n — независимые случайные
величины (СВ), распределенные по Пуассону, равные числу фотонов одного импульса, прошед-
ших без ослабления ОК толщиной P в д.с.п.; Ej, j = 1… N1i(P) — независимые одинаково распре-
деленные СВ, представляющие собой энергии, оставленные в чувствительном слое детектора j-ми
фотонами; Eb i, i = 1…n — независимые одинаково распределенные СВ, являющиеся энергетиче-
скими эквивалентами собственных сигналов детектора; ∆1 — энергетический эквивалент цифро-
вой единицы.
Формула для оценки ∆1 выглядит следующим образом:
k
ADC
∆
1
=C
ADC E
C
(
N
1
(0)E(0)+E
b
) (
2
−1
)
,
(3)
где
N
1
(0)
— среднее количество фотонов, падающих на фронтальную поверхность чувствитель-
ного объема детектора за один импульс без ОК и зарегистрированных в нем; E(0) — средняя
энергия зарегистрированных фотонов без ОК;Eb — среднее значение энергетического эквивален-
та собственных шумов детектора.
Будем исходить из постоянства коэффициента CE на всем диапазоне изменения аналоговых
сигналов. В результате этого допущения и совместного анализа (2) и (3) коэффициент CE может
быть исключен из дальнейшего рассмотрения.
Приведем формулы, уточняющие выражения (2), (3). Для этого используем результаты [18],
приведенные в соответствие с понятием эффективной энергии Eeff [20, 21] применительно к рас-
сматриваемой задаче.
Случайная величина N(P) распределена по Пуассону с параметрами [18]:
Emax
2
N(P)
=σ
N(P)
=
N
f
(
E
,
E
)exp(−
P(E
))
ε
(
E,M
,Z
,p
)
dE
=
1
∫
max
d
d
d
0
(4)
=
N
exp(−P(
E
))ε
E
,
M
,
Z
,
p
,
1
eff
(
eff
d
d
d
)
где P (E ) — толщина ОК в д.с.п. для фотонов с энергией E;
ε
(
E,
M
d
,Z
d
,p
p
)
— эффективность
регистрации детектором с чувствительным элементом с параметрами
M
d
,
Z
d
,p
d
. Формула для
оценки эффективности регистрации имеет вид:
ε
(
E,M
d
,
Z
d
,
p
p
)
=1
−
exp
(
−m
(
E,
M
d
,
Z
d
,
p
d
)
ρ
d d
h
)
,
(5)
здесь
m E,M
d
,Z
d
,p
d
)
— массовый коэффициент ослабления (МКО) фотонов с энергией E:
nd
∑
p
di
M
di
m
(
E
,
Z
di
)
i=1
m
(
E
,
M
,
Z
,
p
)
=
(6)
d
d
d
nd
p
M
∑
di
di
i =1
Дефектоскопия
№ 7
2022
42
В.Ю. Жвырбля, С.П. Осипов, Д.А. Седнев
Случайная величина E(P) характеризуются первым начальным моментом (средним значением
2
E(P)) и вторым начальным моментом
E
(P) [18]. Для полного описания СВ используют плот-
ность распределения вероятностей. Для E(P) плотность распределения g(E(P)) является многопа-
раметрической, то есть зависит от набора параметров (E
max
, P,M
d
,Z
d
,p
d
) и описывается выраже-
нием:
f
(E,E
)exp(−P
(
E))
ε
(
E,M
,Z
,p
)
max
d
d
d
g
(
E(P),E
, P,M
,Z
,p
)
=
max
d
d
d
E
max
(7)
f
(E,E
)exp(−
P
(
E))
ε
(
E
,M
,Z
,p
)
dE
∫
max
d
d
d
0
Формулы для оценки k-го начального момента E(P) выглядят следующим образом:
Emax
k
k
E
(P)
= ∫ E
ab
(E)g
(
E(
P),E
max
, P,
M
d
,
Z
d
,p
d
)
dE,
(8)
0
где Eab(E) — поглощенная энергия зарегистрированного фотона. Эта энергия сложным образом
зависит от материала чувствительного элемента детектора и его размеров [22], поэтому в первом
приближении можно считать, что Eab(E) ≈ E.
Выражения (1)—(8) являются основой для построения имитационной модели формирования
цифровых радиографических изображений для условий ограниченного «квантового голодания» с
целью оценки эффективности подхода, связанного с анализом редких событий. Эти выражения
необходимо дополнить блоком зашумления сигналов в условиях «квантового голодания».
1.3. Моделирование зашумленных сигналов в условиях «квантового голодания»
Из анализа выражения (2), предназначенного для оценки моделирования цифровых сигналов,
можно сделать вывод о применимости прямого моделирования ЦС, разыгрывая для каждого заре-
гистрированного фотона его поглощенную энергию.
Для моделирования ЦС на первом этапе определяется параметр распределения Пуассона
N(P). Распределение Пуассона является дискретным и розыгрыш соответствующей СВ сопряжен
с определенными сложностями [23]. Эти сложности легко устранимы, так как в любых языках
программирования есть встроенные функции для моделирования СВ для наиболее популярных
распределений, в том числе и для распределения Пуассона. Предполагается, что разрабатываемая
математическая модель будет трансформирована в соответствующую имитационную модель с
последующей ее реализацией в системе для математических вычислений MathCad. Широкое при-
менение системы MathCad для имитационного моделирования в физике и приложениях обуслов-
лено близостью языка MathCad к языку математики [11, 18, 24—26]. Для моделирования N(P)
необходимо знать число фотонов от одного импульса N1 [18, 27], попадающих на фронтальную
поверхность чувствительного элемента детектора. С достаточной степенью подробности этот под-
ход изложен в [18]. Точная оценка ожидаемого среднего числа фотонов N1 может быть определена
на основе обработки результатов натурного эксперимента для конкретной модели системы ЦР,
исходя из вычисленной относительной погрешности ЦС для фрагмента откалиброванного радио-
графического изображения свободного от тени ОК.
Для моделирования СВ E(P) вводится плотность распределения (7), после чего эта плотность
преобразуется в функцию распределения и табулируется. Обратная функция к полученной число-
вой зависимости аппроксимируется сплайнами. Данные операции необходимы для розыгрыша
рассматриваемой СВ E(P) с высокой производительностью. Такой способ основан на моделирова-
нии СВ методом обратных функций [28].
Такой подход позволяет моделировать общую СВ — цифровой сигнал D(P).
В процессе моделирования необходимо помнить, что значение параметра Пуассона N (P) не
должно быть меньше единицы. Минимальное же значение цифрового сигнала M [18] устанавли-
вается в процессе вычислительного эксперимента.
В силу отмеченной близости математического языка и языка MathCad приведем алгоритм
имитационного моделирования для условий поставленной задачи максимально близко к тек-
сту MathCad программы. Такой подход позволит избежать недопониманий при изучении тек-
ста работы, и тем более при необходимости повторения результатов вычислительных экспери-
ментов.
Дефектоскопия
№ 7
2022
Повышение проникающей способности систем цифровой радиографии на основе анализа...
43
2. АЛГОРИТМ ИМИТАЦИОННОГО МОДЕЛИРОВАНИЯ ЦИФРОВЫХ
РАДИОГРАФИЧЕСКИХ ИЗОБРАЖЕНИЙ В УСЛОВИЯХ «КВАНТОВОГО ГОЛОДАНИЯ»
За основу алгоритма имитационного моделирования имитационного моделирования цифровых
радиографических изображений в условиях «квантового голодания» были взят ряд результатов и
рекомендаций из [11, 18]. Алгоритм состоит из подготовительного этапа, блока задания базовых
функций, блока описания структуры объекта контроля, подпрограммы формирования цифрового
радиографического изображения с суммированием изображений, подпрограммы визуализации
цифровых радиографических изображений. Здесь будем исследовать позитивные радиографиче-
ские изображения, в которых большей яркости изображения соответствует большая радиационная
прозрачность ОК.
2.1. Подготовительный этап
2.1.1. Формирование таблиц МКО гамма-излучения
Подготовительный этап включает в себя формирование таблиц, содержащих МКО гамма-излу-
чения для энергий в диапазоне от 1 кэВ до 10 МэВ для элементов от водорода (H, Z=1) до урана
(U, Z = 92). Указанный диапазон энергий обусловлен тем, что в индустриальных приложениях
высокоэнергетической цифровой радиографии в качестве источников тормозного излучения при-
менятся бетатроны и линейные ускорители электронов с максимальной энергией Emax до 10 МэВ
[29]. При необходимости диапазон энергий гамма-квантов для таблиц МКО может быть расширен.
Для формирования таблиц МКО использовалась база данных по ослаблению гамма-излуче-
ния [9]. Данная база данных находится в открытом доступе и обладает удобным интерфейсом,
позволяющим пользователю сформировать свой вариант набора данных, учитывающий диапа-
зон энергий гамма-излучения, атомные номера (символы) интересующих химических элемен-
тов, объем информации по видам взаимодействий гамма-излучения с веществом, варианты
выходных единиц (см2/г, барн/атом).
Наиболее корректно база данных [19] работает с броузером Firefox, который позволяет сохра-
нять таблицы МКО в текстовом формате (*.txt). Документы MathCad органично взаимодействуют
с файлами с расширением *.txt. Для хранения таблиц МКО рекомендуется создать специальную
папку, например, X-Ray_Data. Для простоты сопоставления заголовок файла соотносится с атом-
ным номером химического элемента Z, например, файл 1.txt соответствует водороду, а файл 92.txt
— урану. Каждый текстовой файл состоит из заголовка файла и собственно таблицы (двух или
трехмерной матрицы), первый столбец которой — энергии EZ в МэВ, а второй — МКО в см2/г;
третий столбец — массовый коэффициент поглощения в см2/г. Следует отметить, что число строк
в упомянутых матрицах зависит от Z и увеличивается с его ростом, что связано с увеличением
количества характеристических линий в рассматриваемом энергетическом диапазоне.
2.1.2. Ввод и первичная трансформация файлов с МКО
Для удобства использования индексов в MathCad имеется параметр ORIGIN для задания
начального значения индексов, например, ORIGIN:=1, означают, что индексы, нумеруются с еди-
ницы. Здесь «:=» — символ присваивания.
Пусть в папке X-Ray_Data содержатся файлы 1.txt, 2.txt,…, 83.txt и файл ZZM.txt, состоящий из
двух столбцов: первый — атомный номер химического элемента Z; второй — ZM, отношение атом-
ного номера Z к молярной массе M соответствующего химического элемента. Данный файл необ-
ходим для расчета МКО сложных химических веществ. В MathCad имеются средства для работы
со строковыми переменными, что облегчает формирование списков Ввода/Вывода в соответству-
ющих программах. Строковые переменные в MathCad заключатся в кавычки — ″Строковая пере-
менная″. Применительно к рассматриваемому случаю процесс ввода всех необходимых файлов
описывается следующей конструкцией:
Z
:=
n m
file
:=
concat
(
"X -Ray_Data
/",num2str(Z),".txt"
)
Z
,
(9)
C
:=
READPRN
(
file
)
Z
Z
ZM
:=
READPRN("X -Ray_Data
/
ZZM txt
")
Дефектоскопия
№ 7
2022
44
В.Ю. Жвырбля, С.П. Осипов, Д.А. Седнев
здесь fileZ — имя файла для вещества с атомным номером Z; «n..m» — изменение переменной в
диапазоне от n до m; concat — функция объединения строковых переменных, перечисленных через
запятую; num2str(Z) — преобразование числа Z в строковую переменную.
Системы ЦР и РКТ используются в разных энергетических диапазонах — от десятков кэВ до
10 МэВ, поэтому шкала энергий EZ переводится в кэВ. Первичная трансформация сводится допол-
нительно к введению МКО, атомных номеров элементов Ze и молярных масс Mz:
1
2
1
2
E
:= 1000
(
C
)
m
:=
(
C
)
Ze:=
(
ZM
)
Mz:=
(
ZM
)
,
(10)
Z
Z
Z
Z
здесь j
— j-й столбец матрицы.
В ряде случаев необходимо расположить энергии по возрастанию, а в таблицах [19] присут-
ствуют дублированные характеристические линии Eχ, для которых m(Eχ-) ≠ m(Eχ+), поэтому необ-
ходима корректировка вектора
E
в упорядоченный вектор EcZ:
Z
γ ← 0.00001
i0 ←
rows
(
E
)
Z
Ec
Z
:=
for i∈2..i0
,
(11)
(
E
)
←
(
E
) ⋅(
1+γ
)
if
(
E
)
=
(
E
)
Z
i
Z
i
Z
i
Z
i−1
EZ
здесь γ — параметр, меньший относительной погрешности оценки энергии характеристического
излучения; |— символ подпрограммы; ← — символ присваивания в подпрограммах; rows — коли-
чество строк в матрице; for — оператор цикла; i ∈ 2..i0 — диапазон изменения переменной i от 2
до i0.
2.1.3. Интерполяция зависимостей МКО гамма-излучения от энергии
Система MathCad обладает широкими возможностями для интерполяции, но в рассматривае-
мом случае достаточно линейной интерполяции (встроенная функция linterp). Формула для интер-
поляции МКО от энергии гамма-кванта x для элемента с атомным номером Z имеет вид:
mas x,
Z):= linterp
(
E
,
m
, x
)
(12)
Z
Z
Выражение (12) может быть использовано для интерполяции МКО для чувствительного объ-
ема детектора masd от энергии x. В соответствии с (6) интерполяция masd(x) применительно к
системе MathCad запишется следующим образом:
nd
nd
masd x):
pd
Md
m x,Zd
)
pd
Md
(13)
=∑
i
i
i
∑
i
i
i
=1
i
=1
Векторы Md, Zd, pd в (13) совпадают с соответствующими векторами из (6)
M
,Z
,p
d
d
d
M
Ze
1
53
53
Для сцинтиллятора CsI —
Md
:=
,
Zd
:=
,
pd
:
=
M
Ze
1
55
55
2.2. Блок задания базовых функций
2.2.1. Эффективность регистрации
Функция для вычисления зависимости эффективности регистрации ε от энергии x, реализую-
щая (5) в программе MathCad, имеет вид:
ε x):=
1
−
exp(
−masd x)
ρ
h
)
(14)
d d
Дефектоскопия
№ 7
2022
Повышение проникающей способности систем цифровой радиографии на основе анализа...
45
2.2.2. Числовой энергетический спектр
Для высокоэнергетического тормозного излучения распределение интенсивности I по энергии
фотонов x — dI/dx близко к прямой линии [30, 31]:
2
Shiff
(
Em, x
)
:= 2
(
Em - x
)
Em
,
(15)
здесь Em — максимальная энергия тормозного излучения; x — энергия фотонов.
Для последующего моделирования необходим числовой энергетический спектр f(Em, x). Так
как dI(x)/dx = x · dN(x)/dx = Shiff(Em, x), то соответствующая функция имеет вид:
Shiff
(
Em,
x
)
f
(
Em x):=
,
(16)
Em
Shiff
(
Em
, x
)
x
dx
∫
x
E
0
здесь E0=1 кэВ — минимальное значение энергии фотонов в таблицах МКО из [19].
2.2.3. Зависимость потока частиц и интенсивности от параметров ОК
На практике низкий уровень радиационной прозрачности проявляется при измерениях прони-
кающей способности Hlim. Эти измерения проводят на стальных тестовых образцах с плотностью
ρst=7,56 г/см3, поэтому примем:
mst(x) := mas(x,26), ρst := 7,56.
(17)
Ожидаемое число фотонов N0 на фронтальной поверхности детектора при отсутствии ОК
вычисляется с учетом апертуры детектора ad×bd мм2, количества интегрируемых аналоговым
образом импульсов kimp, числа фотонов NN0, падающих на 1 см2 поверхности детектора.
Детекторы ИДК ТПУ [31] характеризуются параметрами:
5
ad := 6, bd := 5, kimp := 1, NN 0 := 2,4 ⋅10
(18)
С учетом сказанного выражение для оценки N0 выглядит следующим образом:
N0:= NN0⋅ad ⋅bd ⋅kimp⋅0,01.
(19)
Поток фотонов N (Em, Hst,hd ) вычисляется по формуле:
Em
−mst(x)⋅ρst⋅Hst
N(Em,Hst,hd):=
floorN0⋅
f
(
Em x)⋅e
⋅ε(x,hd)dx
,
(20)
∫
E
0
где floor(q) — функция MathCad, целая часть числа q.
Средняя энергия Emean(Em, Hst,hd ) зарегистрированного фотона оценивается следующим
образом:
Em
−moc x)⋅ρstHst
x⋅ f
(
Em x)
⋅ e
⋅ε(x,hd)dx
∫
E0
Emean(Em,Hst,
hd):
=
(21)
Em
−moc(x
)⋅ρstHst
f
(
Em x)
⋅
e
⋅ε(x,hd)dx
∫
E
0
В свою очередь интенсивность энергии I (Em, Hst,hd ) оценивается по формуле:
I(Em,Hst,hd):= Emean(Em,Hst,hd)⋅N(Em,Hst,hd).
(22)
В соответствии с (5) выражение для расчета аналога цифровой единицы имеет вид:
kadc
∆(Em,hd,
kadc):=I
(
Em
,0,hd
)
(
Cadc ⋅
(
2
−1
))
,
(23)
здесь kadc, Cadc — разрядность АЦП и коэффициент заполнения диапазона изменения ЦС.
Дефектоскопия
№ 7
2022
46
В.Ю. Жвырбля, С.П. Осипов, Д.А. Седнев
2.3. Описание объектов контроля
2.3.1. Общие положения
При описании ОК воспользуемся подходом, предлагаемым в работах [11, 32]. Данный подход
основан на представлении ОК в качестве объединения непересекающихся фрагментов, причем в
ряде случаев эти фрагменты имеют правильную форму и ориентацию. Отмеченный факт характе-
рен для некоторых тестовых объектов (ТО) (эталонов), предназначенных оценки параметров
систем ЦР. К таким ТО относятся [11, 32], например, эталоны для оценки проникающей способ-
ности, эталоны чувствительности, канавочные эталоны для оценки пространственного разреше-
ния, ступенчатые эталоны и т.п. Все эти эталоны изготавливают, как правило, из стали. Для подоб-
ных ТО ограниченность количеств фрагментов в сочетании с их формой (прямоугольные паралле-
лепипеды) приводит и к ограниченному набору толщин по стали.
Эффект «квантового голодания» проявляется для ОК значительной толщины, поэтому в каче-
стве примера можно ограничиться ступенчатым объектом. Такой объект может быть сформирован
стальными пластинами из эталона для оценки проникающей способности. Пусть толщина пластин
Hr кратна 6 мм, пластины имеют поперечные размеры Ar×Ar = 600×600 мм2. Ступенчатый объект
формируется параллельным смещением пластин друг относительно друга. Для описания кон-
струкции ОК введем систему координат XOY, где координату x свяжем с направлением сканирова-
ния, а координату y с линейкой детекторов. Пусть границы пластин параллельны осям координат,
а центры пластин расположены по линии y = y0 со смещением ∆x.
2.3.2. Функции и подпрограммы
Для удобства применения подхода к описанию ОК [32] введем логическую функцию W(x, y,
x0, y0) принадлежности точки (x, y) пластине с центром в точке (x0, y0):
W(x, y,x0, y0):=
(
x-x0
≤ Ar 2 ∧(
y-y0
≤ Ar 2
)
(24)
Подпрограмма для оценки толщины ОК с учетом окружающего воздуха над плоскостью фор-
мируемого изображения HST(x, y, x0, y0) имеет следующий вид:
kr0←
rows(x0)
for kr
∈1..kr0
Hrr
kr
←
Hr if W x,y,x0 ,
kr
y0)
H x,y,x0,y0):=
,
(25)
Hrr
←0
otherwise
kr
kr
0
∑Hrrkr
kr
=1
здесь otherwise — оператор MathCad «в противном случае» дополняет пространство событий в
операторе if.
Воспользуемся ограниченностью количества уровней толщин и введем соответствующую
функцию LevelH(x, y, x0, y0):
L
←
1
for kr∈1..kr0
LevelH x,y,x0,y0):=
(26)
L
←
kr
+1
if W x,y,x0 ,
y0)
kr
L
Конструкция (26) является исключительно простой, но при вычислительном моделировании
цифровых радиографических изображений крупногабаритных объектов она обеспечивает ускоре-
ние процесса вычислений в сотни раз.
Каждому уровню толщины LH соответствует свои уровни числа фотонов LN, среднего значе-
ния энергии LE и аналогового сигнала LI. Отметим, что общее количество уровней на единицу
больше количества ступеней kL0:
Дефектоскопия
№ 7
2022
Повышение проникающей способности систем цифровой радиографии на основе анализа...
47
LH
:=
(
kL
−1)Hr,
LN
:=
N(Em,LH
,hd ),
kL
kL
kL
(27)
LE
:=
Emean(Em,LH
,hd
),
LI
:=
LN
⋅
LE
kL
kL
kL
kL
kL
Введенных выше функций, подпрограмм и векторов достаточно для перехода к подпрограмме
формирования цифровых радиографических изображений.
2.3. Подпрограмма формирования цифрового радиографического изображения
Входными параметрами подпрограммы XrayImage являются поперечные размеры ОК Xoc×Yoc,
поперечные размеры детектора ad×bd, количество суммируемых изображений n. Заранее должны
быть определены среднее значение mb и среднеквадратическое отклонение σb темнового сигнала
детектора.
Подпрограмму приведем в упрощенной версии без розыгрыша энергий зарегистрированных
фотонов. На момент вызова все параметры подпрограммы XrayImage должны быть определены.
Xoc
floor
NX
ad
←
NY
Yoc
floor
bd
for ix∈1..NX
+1
for iy∈
1..NY
+1
Xoc
−
+
(
ix
−1
)
⋅ ad
X
ix
2
←
Y
Yoc
iy
−
+
(
iy
−1
)
⋅bd
2
lh
←
LevelH
(
x y,x0,
y0
)
Noc← LNlh
Eoc← LElh
rpois
(
n,
Noc
)
if Noc
≥1
XrayImage(Xoc,Yoc,ad,
bd n):=
r
←
,
(28)
null n)
otherwise
black ← rnorm
(
n
,
mb
,
σb
)
n
Eoc⋅r
i
ID
←
floor
+black
−mb
iy
,ix
∑
i
i=1
∆
1
for iy∈1..NY
+1
iw
y ix
i=1
MD
←
iy
iw
for ix∈1..
NX
+1
for iy
∈
1..NY
+1
max(0.5,
ID
)
iy
,ix
DH
←-
ln
iy
,ix
MD
iy
ID
Дефектоскопия
№ 7
2022
48
В.Ю. Жвырбля, С.П. Осипов, Д.А. Седнев
здесь rpois(n, Noc), rnorm(n, mb, σb)— программы розыгрыша векторов СВ, распределенных по
Пуассону и Гауссу, в данном случае векторы имеют размерность n; null(n) — нулевой вектор раз-
мерности n.
2.4. Подпрограмма формирования полутоновых изображений
На вход подпрограммы визуализации VizD поступают: цифровое радиографическое изображе-
ние D := XrayImage( Xoc,Yoc,ad,bd,n); минимальный и максимальный уровни тона Pmi, Pma,
удобные для восприятия; вектор pp — цифровых координат в системе RGB для раскраски фоновой
части изображения; вектор pr — цифровых координат в системе RGB для раскраски части
изображения с низким уровнем ЦС. В качестве векторов pp, pr могут быть использованы
1
1
pp:=
1 ,
pr:=
0
. В результате фон будет изображен в оттенках серого, а пиксели с низким
1
0
уровнем ЦС будут выделены красным цветом.
Nx
cols(D)
←
Ny
rows(D)
Dmi
min(D
)
←
Dma
max(D
)
for ix
∈1..NX
for iy∈1..NY
D
iy,ix
-
Dmi
VizD(D,sh,Pmi,Pma,
pp,
pr):=
P
←
Pma- floor
(
Pma−Pmi
)
,
(29)
iy,ix
Dma-Dmi
riy,
ix
kadc
-D
iy,ix
g
←
pp⋅P
if Cadc⋅
2
−1
⋅e
>1
iy
,
ix
iy,ix
(
)
b
iy
,ix
riy,
ix
g
←
pr
⋅ 255
otherwise
iy
,
ix
biy,
ix
augment r,g,b)
здесь min(D), max(D) — минимальное и максимальное значения массива D, которые находятся с
помощью встроенных функций MathCad; augment — функция объединения матриц (векторов) в
общую матрицу.
2.5. Вывод и визуализация полутоновых изображений
В системе MathCad для ввода/вывода изображений в оттенках серого используются функции
READBMP/WRITEBMP, а в оттенках цветного — READRGB/WRITERGB. Для рассматриваемой
задачи более предпочтительной является вторая группа функций.
На выходе подпрограммы VizD формируется не само изображение, а матрица VI, ассоциируе-
мая с ним VI := VizD(D,sh, Pmi, Pma, pp, pr) . Для превращения матрицы VI в цветное изображение
вначале необходимо сформировать текстовую переменную fileI, содержащую информацию о месте
записи файла (изображения), его формате, количестве суммируемых цифровых изображений, а
затем выполнить соответствующую операцию ввода. Упомянутый комплекс действий описывается
двумя командами:
Дефектоскопия
№ 7
2022
Повышение проникающей способности систем цифровой радиографии на основе анализа...
49
fileI
:=
concat
("
BMP
/","
Sum
",
nim2
str n),"
digital image
",".bmp
")
(30)
(
WRITERGB
fileI
):=
VI
Папка BMP предназначена для хранения изображений и создается заранее.
Приведенный алгоритм (1)—(30) легко трансформируется в программу в системе MathCad с
уточнением входных данных.
3. РЕЗУЛЬТАТЫ МОДЕЛИРОВАНИЯ
3.1. Исходные данные
Часть данных уже демонстрировалась выше в некоторых блоках и этапах разработанного алго-
ритма. Перечислим полный набор исходных данных.
Источник высокоэнергетического тормозного излучения — малогабаритный бетатрон ИДК
ТПУ [31] с максимальной энергией Emax = 9 МэВ. Оценочное и подтвержденное экспериментально
количество фотонов на 1 см2 на расстоянии 4,2 м от источника, приведенное к одному импульсу,
N1 = 240000 (30 сГр/мин).
Линейка детекторов: материал сцинтилляторов CdWO4; поперечные размеры детекторов
ad×bd
= 6×5 мм2; толщина детекторов hd = 35 мм; разрядность АЦП kADC = 16 бит; коэффициент
сужения диапазона цифровых сигналов CADC = 0,8. Коэффициент геометрического увеличения
равен 3.
Базовый объект контроля: стальной ступенчатый объект, набранный из пластин размером
600×600 мм2; ширина объекта по оси y 600 мм; высота ступени 24 мм; количество ступеней 13;
ширина ступени 30 мм. Число ступеней базового ОК было выбрано из условия, что число фотонов,
регистрируемых детектором, не должно быть меньше 1. Это ограничение связано с рассматривае-
мой концепцией моделирования. Для условий моделирования максимальная толщина объекта по
стали 312 мм.
Выше отмечено, что существенное несоответствие диапазонов изменения АС и ЦС может при-
вести к эффекту, аналогичному эффекту «квантового голодания», несмотря на сравнительно высо-
кий уровень числа регистрируемых фотонов. Этот эффект напоминает по признакам «квантовое
голодание», назовем его эффектом «псевдо-квантового голодания». Отмеченный эффект характе-
ризуется упомянутым выше параметром CADC. В соответствии со сказанным выше уменьшим
высоту ступени в базовом ОК до 18 мм, что приведет к уменьшению максимальной массовой
толщины ОК до 234 мм по стали.
3.2. Результаты моделирования цифровых радиографических изображений
3.2.1. Базовый ступенчатый объект, условия ограниченного «квантового голодания»
На рис. 1 приведены цифровые радиографические изображения испытуемого ступенчатого
объекта с вариацией числа исходных изображений от 1 до 1024.
Сравнение изображений, приведенных на рис. 1, подтверждает эффективность цифрового сум-
мирования радиографических изображений в условиях ограниченного «квантового голодания».
3.2.2. Базовый ступенчатый объект, условия «псевдоквантового голодания»
Для исследования эффекта «псевдоквантового голодания» выберем ступенчатый объект, опи-
санный выше, с существенно меньшим значением параметра CADC.
На рис. 2 приведены цифровые радиографические изображения исследуемого ОК с максималь-
ной толщиной 234 мм по стали для n = 1; 4; 16 и значения параметра CADC = 0,1.
В результате сравнения изображений, представленных на рис. 1 и 2, можно подтвердить схо-
жесть эффектов «квантового голодания» и «псевдоквантового голодания». Отмечается повышение
качества изображений с увеличением параметра n.
Следует отметить, что эффект «псевдоквантового голодания» обуславливается также недоста-
точно высокой разрядностью АЦП (см. (3)).
Дефектоскопия
№ 7
2022
50
В.Ю. Жвырбля, С.П. Осипов, Д.А. Седнев
Рис. 1. Цифровые радиографические изображения ступенчатого объекта с вариацией числа суммируемых изображений n.
Эффект «квантового голодания».
Рис. 2. Цифровые радиографические изображения ступенчатого объекта с вариацией числа суммируемых изображений n.
Эффект «псевдоквантового голодания».
3.2.3. Эталон проникающей способности
При формировании изображений эталонов для оценки проникающей способности систем ЦР
целенаправленно создается эффект «квантового голодания» [33].
При моделировании изображений эталонов для оценки проникающей способности были
использованы формулы из [34] для корректировки пространственного описания ОК.
На рис. 3 приведены изображения набора прямоугольных стальных пластин толщиной
Hst = 324 мм с «наконечником стрелы» из стали толщиной Ha = 84 мм для двух положений наконеч-
ника, Hlim= Hst+Ha = 408 мм.
Анализ данных, приведенных на рис. 3, подтверждает вывод об эффективности цифрового
суммирования изображений в условиях «квантового голодания». Наблюдается улучшение четко-
сти границ «стрелы наконечника» с ростом параметра n.
В работе [27] отмечена возможность достижения проникающей способности Hlim = 450 мм по стали
для тормозного излучения с Emax = 9 МэВ при условии попадания не менее 5 фотонов на пиксель раз-
мером 10×10 мм2 на расстоянии 7,5 м. Из сопоставления параметров систем ЦР из [27] и [31] следует,
что значение параметра N1 для упомянутой в условии системы [27] примерно в 26 раз больше, чем для
Дефектоскопия
№ 7
2022
Повышение проникающей способности систем цифровой радиографии на основе анализа...
51
Рис. 3. Цифровые радиографические изображения эталона проникающей способности с вариацией числа
суммируемых изображений n. Эффекты «квантового голодания» — Hlim = 408 мм и «псевдоквантового голодания»
для условий [27] — Hlim=450 мм.
той, что использовалась в модельном примере. Для такого уровня N1 был сформированы изображения,
приведенные на рис. 3. Отмеченное в [27] условие позволяет объяснить наблюдаемый эффект «псевдо-
квантовым голоданием» — не полным соответствием диапазонов изменения АС и ЦС.
ЗАКЛЮЧЕНИЕ
Разработана математическая модель формирования высокоэнергетических цифровых радио-
графических изображений в условиях низких уровней цифровых сигналов. Модель предназначена
для повышения эффективности анализируемых систем применительно к контролю крупногаба-
ритных объектов за счет цифрового суммирования изображений. На основе модели разработан
алгоритм и программа на MathCad, предназначенные для симуляции цифровых изображений
испытуемых объектов. На основе анализа результатов вычислительных и натурных экспериментов
проиллюстрирована принципиальная возможность повышения эффективности систем высокоэ-
нергетической цифровой радиографии на основе учета редких событий.
Работа выполнена при финансовой поддержке Министерства науки и высшего образования
Российской Федерации в рамках Государственного задания «Наука», проект № FSWW-2020-0014.
СПИСОК ЛИТЕРАТУРЫ
1. Du Plessis A., MacDonald E., Waller J.M., Berto F. Non-destructive testing of parts produced by laser
powder bed fusion // Fundamentals of Laser Powder Bed Fusion of Metals. Elsevier, 2021. P. 277—300.
2. Bossema F.G., Domínguez-Delmás M., Palenstijn W.J., Kostenko A., Dorscheid J., Coban S.B., Her-
mens E., Batenburg K.J. A novel method for dendrochronology of large historical wooden objects using line
s41598-021-90135-4
Дефектоскопия
№ 7
2022
52
В.Ю. Жвырбля, С.П. Осипов, Д.А. Седнев
3. Nyberg J.V.E., MacGregor D., Azaiez F. 5. Complementary methods: γ-beam techniques, X-ray
fluorescence (XRF) and nuclear magnetic resonance (NMR) // Nuclear physics for cultural heritage. EDP
Sciences. 2021. P. 37—53.
4. Karikh V.P., Pevchenko B.V., Kurbatov A.V., Okhotnikov A.A., Skokov A.A. Three-dimensional X-ray
tomography of industrial objects with limited sizes of the recording screen // Russian Journal of
[Карих В.П., Певченко Б.В., Курбатов А.В., Охотников А.А., Скоков А.А. Трехмерная рентгеновская
томография промышленных объектов при ограниченных размерах регистрирующего экрана //
Дефектоскопия. 2021. № 9. С. 54—62.]
5. Ewert U. Current developments in digital radiography and computed tomography from nm to macro
scale / Proceedings of the 12th European Conference on Non-Destructive Testing (ECNDT 2018), Gothenburg,
19.04.2022).
6. Movafeghi A., Mohammadzadeh N., Yahaghi E., Nekouei J., Rostami P., Moradi G. Defect detection of
industrial radiography images of ammonia pipes by a sparse coding model // Journal of Nondestructive Evaluation.
7. Kaczmarek Ł., Kozłowska A., Maksimczuk M., Wejrzanowski T. The use of X-ray computed
microtomography for graptolite detection in rock based on core internal structure visualization // Acta Geologica
8. Udod V.A., Osipov S.P., Wang Y. Estimating the influence of quantum noises on the quality of material
identification by the dual-energy method // Russian Journal of Nondestructive Testing. 2018. V. 54. No. 8.
квантовых шумов на качество распознавания материалов методом дуальных энергий // Дефектоскопия.
2018. № 8. С. 50—65.]
9. Li L., Li R., Zhang S., Zhao T., Chen Z. A dynamic material discrimination algorithm for dual MV energy
X-ray digital radiography // Applied Radiation and Isotopes. 2016. V. 114. P. 188—195. https://doi.
org/10.1016/j.apradiso.2016.05.018
10. Hall C. Large-scale imaging on the Australian Synchrotron imaging and medical beamline (IMBL) //
11. Osipov S.P., Chakhlov S.V., Kairalapov D.U., Sirot’yan E.V. Numerical modeling of radiographic
images as the basis for correctly designing digital radiography systems of large-sized objects // Russian Journal
[Осипов С.П., Чахлов С.В., Кайралапов Д.У., Сиротьян Е.В. Численное моделирование радиографиче-
ских изображений — основа корректного проектирования систем цифровой радиографии крупногаба-
ритных объектов // Дефектоскопия. 2019. № 2. С. 43—55.]
12. Busi M., Kehl C., Frisvad J.R., Olsen U.L. Metal artifact reduction in spectral X-ray CT using spectral
13. Zeng G.L. Photon Starvation Artifact Reduction by Shift-Variant Processing // IEEE Access. 2022.
14. Morin R.L., Seibert J.A. Considerations for selecting a digital radiography system // Journal of the
15. Zolkin A.L., Munister V.D., Bogaevskaya O.Y., Yumashev A.V., Kornetov A.N. Promising directions of
information technologies applications for improving of control and measuring devices in medical industry //
Journal of Physics: Conference Series. — IOP Publishing, 2021. V. 2094. No. 3. No. article 032001. https://
doi.org/10.1088/1742-6596/2094/3/032001
16. Andriole K.P. Medical imaging modalities and digital images / Practical Imaging Informatics. Springer,
17. Beckmann M., Bhandari A., Krahmer F. The Modulo Radon Transform: Theory, Algorithms, and
Applications // SIAM Journal on Imaging Sciences. 2022. V. 15. No. 2. P. 455—490. https://doi.
org/10.1137/21M1424615
18. Osipov S.P., Usachev E.J., Chakhlov S.V., Schetinkin S.A., Osipov O.S. Inspection of bulk cargoes and
liquids by the dual energy method // Radiation Physics and Chemistry. 2020. V. 177. No. article 109133.
19. Berger M.J., Hubbell J.H., Seltzer S.M., Chang J., Coursey J.S., Sukumar R., Zucker D.S., Olsen K.
XCOM: Photon cross sections database, NIST standard reference database 8 (XGAM). https://doi.
org/10.18434/T48G6X
20. Petrushanskii M.G., Kornev E.A. On determination of the effective energy of a mixed X-ray beam //
[Петрушанский М.Г., Корнев Е.А. К определению эффективной энергии смешанного пучка рентгенов-
ского излучения // Медицинская техника. 2006. № 2. С. 46—49.]
21. Zhang D., Lai W., Wu J., Du H., Zhao R., Fan S. Establishment of ISO 4037-1 X-ray narrow-spectrum
joe.2018.9126
Дефектоскопия
№ 7
2022
Повышение проникающей способности систем цифровой радиографии на основе анализа...
53
22. Zav’yalkin F.M., Osipov S.P. Dependence of the mean value and fluctuations of the absorbed energy on the
23. Barbiero A., Ferrari P.A. Simulation of correlated Poisson variables // Applied Stochastic Models in
24. Zhao N., Chen X. Exploring the integration of MathCad-assisted mathematics experiments into higher
mathematics teaching // Computer-Aided Design and Applications. 2022. 19(S1). P. 117—127. https://doi.
org/10.14733/cadaps.2022.S1.117-127
25. Sarantites D.G., Sobotka L.G. Computational templates for introductory nuclear science using
26. Chang J., Soldatov A., Sorokin P., Kostina M., Soldatov A., Xu S., Chiriev A. Gaussian filtered small-
angle ultrasonic computed tomography using linear arrays / Progress in Material Science and Engineering.
27. Kutsaev S., Agustsson R., Arodzero A., Boucher S., Hartzell J., Murokh A., O’Shea F., Smirnov A.Y.
Electron accelerators for novel cargo inspection methods // Physics Procedia. 2017. V. 90. P. 115—125. https://
doi.org/10.1016/j.phpro.2017.09.036
28. Graham C., Talay D. Strong law of large numbers and Monte Carlo methods / Stochastic Simulation
39363-1_2
29. Kutsaev S.V. Advanced technologies for applied particle accelerators and examples of their use //
ев С.В. Перспективные технологии для прикладных резонансных ускорителей заряженных частиц и
примеры их использования (Обзор) // Журнал технической физики. 2021. Т. 91. № 2. С. 173—206.]
30. Schiff L.I. Nonlinear meson theory of nuclear forces. II. Nonlinearity in the meson-nucleon coupling
31. Khasanov S., Suleymanov R., Safarov A., Safarov A., Tedila H.M., Muratov R. Calculation of exposure
and absorbed dose at irradiation samples by electron bremsstrahlung // Radiation Physics and Chemistry. 2021.
laboratory/rknl/eng/products/iDK (Обращение 27.04.2022).
33. ANSI N42.46. American national standard for determination of the imaging performance of X-ray and
gamma-ray systems for cargo and vehicle security screening, 2008.
34. Osipov S.P., Yadrenkin I.G., Chakhlov S.V., Osipov O.S., Usachev E.Y. Simulation modeling in digital
radiography with allowance for spatial outlines of test objects // Russian Journal of Nondestructive Testing.
Чахлов С.В., Осипов О.С., Усачев Е.Ю. Имитационное моделирование в цифровой радиографии с уче-
том пространственных форм объектов контроля // Дефектоскопия. 2020. № 8. С. 35—48.]
Дефектоскопия
№ 7
2022
