УДК 620.179.152
ВОССТАНОВЛЕНИЕ СИГНАЛОВ СО СЦИНТИЛЛЯЦИОННЫХ ДЕТЕКТОРОВ
© 2022 г. С.П. Осипов1, 2, *, С.А. Щетинкин2, Е.Ю. Усачёв2, С.В. Чахлов1, 2, **, О.С. Осипов3
1Национальный исследовательский Томский политехнический университет,
Россия 634028 Томск, пр-т Ленина, 30
2ООО Диагностика-М, Россия 109316 Москва, Волгоградский пр-т, 42
3Medialooks, Россия 236016 Калининград, ул. Александра Невского, 59
E-mail: *osip1809@rambler.ru; **chakhlov@tpu.ru
Поступила в редакцию 13.04.2022; после доработки 22.04.2022
Принята к публикации 27.04.2022
Предложена математическая модель искажения исходного радиометрического сигнала сцинтилляционным
детектором для случая регистрации непрерывного рентгеновского излучения в процессе формирования цифровых
радиографических изображений. Модель учитывает распределение исходного сигнала по времени, время регистра-
ции излучения, разрядность АЦП, временные и амплитудные параметры быстрой и медленной составляющих про-
цесса высвечивания в сцинтилляционном детекторе. Разработан алгоритм восстановления сигнала после сцинтил-
лятора, основанный на предварительной информации о параметрах сцинтилляционного детектора для непрерывно-
го рентгеновского излучения. Предложены алгоритмы оценки амплитудных и временных параметров быстрой и
медленной составляющих процесса высвечивания в сцинтилляционном детекторе для источников рентгеновского
излучения. Проведена серия вычислительных экспериментов и на основе анализа ее результатов продемонстриро-
вана эффективность разработанных алгоритмов.
Ключевые слова: рентгеновское излучение, цифровая радиография, сцинтилляционные детекторы, быстрые и
медленные компоненты высвечивания, математическое моделирование, имитационное моделирование, восстановле-
ние сигналов после сцинтиллятора.
DOI: 10.31857/S0130308222060057; EDN: BNAJHH
ВВЕДЕНИЕ
Необходимость повышения качества контроля объектов различной природы методами циф-
ровой радиографии (ЦР) и рентгеновской компьютерной томографии (РКТ) обуславливает
настоятельную потребность в совершенствовании способов формирования соответствующей
первичной информации с учетом особенностей процессов регистрации рентгеновского излуче-
ния. Сцинтилляторы широко применяются в качестве первичных преобразователей радиации в
световой поток с последующей его регистрацией фотоприемником и трансформацией получаю-
щихся аналоговых сигналов в цифровые [1—3]. Преобразование энергии рентгеновского излуче-
ния, поглощенного в сцинтилляторе, в энергию светового потока (высвечивание) является инер-
ционным [1, 3]. Инерционность высвечивания приводит к искажению исходных сигналов сцин-
тилляционным детектором и, следовательно, к ухудшению пространственного разрешения
систем ЦР и РКТ [4—6]. Учет этого фактора в алгоритмах формирования обработки информа-
ции в системах ЦР и РКТ, а также в системах досмотрового контроля с функцией распознавания
материалов методом дуальных энергий невозможен без полной информации о параметрах
быстрой и медленной составляющих процесса высвечивания энергии в сцинтилляторе [1, 7—9].
Временное изменение сигнала на выходе фотопреобразователя описывают экспонентой либо
суммой экспонент [10—14]. В научной литературе не в полной мере обсуждены вопросы, свя-
занные с определением параметров быстрой и медленной составляющих процесса высвечива-
ния и восстановлением сигналов сцинтилляционным детектором для источников непрерывного
рентгеновского излучения.
1. МОДЕЛЬ ИСКАЖЕНИЯ СИГНАЛА СЦИНТИЛЛЯЦИОННЫМ ДЕТЕКТОРОМ
1.1. Представление исходного сигнала
В качестве источников рентгеновского излучения в системах ЦР и РКТ используют, как прави-
ло, рентгеновские аппараты непрерывного действия [15].
Остановимся на формировании исходных радиографических изображений в системах сканиру-
ющего типа с линейным регистратором рентгеновского излучения. Пусть ось ОХ связана с направ-
Восстановление сигналов со сцинтилляционных детекторов
49
лением сканирования, а ось OY — с линейкой детекторов. Объект контроля перемещается относи-
тельно системы, состоящей из источника рентгеновского излучения и линейного детектора, с
постоянной скоростью V.
Будем считать сигнал I со сцинтилляционного детектора «идеальным», если он зависит
исключительно от ослабления излучения объектом контроля. Для интенсивных источников
рентгеновского излучения непрерывного действия традиционным режимом регистрации являет-
ся интегральный, поэтому I можно интерпретировать как количество энергии фотонов, передан-
ной детектору.
Будем придерживаться модели формирования цифрового радиографического изображения [16,
17]. Согласно этой модели, сигнал I в момент времени t определяется текущей толщиной ОК в
длинах свободного пробега P:
I(t) = I
0
(t)exp(-P(t)),
(1)
где I0(t) — зависимость энергии рентгеновского излучения, переданной детектору, от времени t без
объекта контроля. Для непрерывного режима излучения соблюдается равенство
0
I t)
C
=I
1.2. Искажение исходного сигнала сцинтилляционным детектором
Пусть в момент времени p в результате воздействия на сцинтиллятор ультракороткого импуль-
са фотонов высвечиваются световые фотоны с общей энергией E, E = CE I. Согласно [10—14],
процесс высвечивания описывается функцией h(t), равной сумме n экспоненциальных функций.
Для n = 2 функция h(t) представима в виде:
0,
t
<
p
h(t)
=
,
(2)
A
exp
( (
t
−
p
)
τ
)
+
A
exp
( (
t
−
p
)
τ
)
,
t
≥
p
f
f
s
s
где Af, As — амплитудные параметры быстрой и медленной составляющих; τf, τs — временные
параметры быстрой и медленной составляющих.
Будем считать, что суммарная энергия светового излучения E высвечивается без потерь.
Отсюда следует, что
∞
A
exp
(
−(
t
−
p
)
τ
)
+
A
exp
(
−(
t
−
p
)
τ
dp= A
τ
+
A
τ
(3)
∫
(
f
f
s
s
))
f
f
s s
p
Введем обозначение:
a
=A
(
A
τ
+A
τ
),
a
=A
(
A
τ
+A
τ
).
(4)
f
f
f
f
s s
s
s
f
f
s s
С учетом (2) и (4) искажение J(t) исходного сигнала I(t) сцинтилляционным детектором описы-
вается выражением:
t
t
J t)
=C
E
I p)h(
t-p)d
p=C
E
I p)
(
a
f
exp
(
−
(t-p)
τ
f
)
+a
s
exp
(
−
(t-p)
τ
s
))
d
p
(5)
∫
∫
0
0
Полученные выражения носят общий характер и могут быть использованы для источников
рентгеновского излучения непрерывного действия.
Совокупность выражений (1)—(5) является математической моделью трансформации (искаже-
ния) идеальных сигналов для источников рентгеновского излучения сцинтилляционным детекто-
ром, учитывающей быструю и медленную составляющие высвечивания.
Данная математическая модель может быть преобразована в соответствующую имитационную
модель с последующей ее реализацией, например, в системе для математических вычислений
MathCad.
Дефектоскопия
№ 6
2022
50
С.П. Осипов, С.А. Щетинкин, Е.Ю. Усачёв и др.
1.3. Интегрирование, дискретизация и оцифровка сигнала на выходе сцинтилляционного
детектора
Аналоговый сигнал на выходе сцинтилляционного детектора J(t) интегрируется за время ∆t,
поступает на вход аналого-цифрового преобразователя (АЦП) с разрядностью kADC, где и оцифро-
вывается.
Пусть аналоговый сигнал J(t) определен на интервале от 0 до tmax. Допустим, что интегратор
обладает свойством мгновенного сброса сигнала, тогда количество точек в дискретизации сигнала
J(t) равно jmax=[tmax/∆t], здесь и далее [x] — целая часть x. Дискретизация сигнала J(t) будет пред-
ставлять собой вектор аналоговых величин
Jd =
(
Jd
, Jd
, Jd
,..., Jd
)
:
1
2
3
j
max
j∆t
Jd
=
J t)dt,
j
=1...j
(6)
j
∫
max
(
j
−1)
∆t
Аналоговые величины
Jd
, Jd
, ... ,
Jd
преобразуются в цифровые
D
, D
,..., D
:
1
2
j
max
1
2
j
max
C
Jd
ADC j
max(Jd)
D
=
,
j
=1 ...
j
,
∆
=
,
(7)
j
max
1
k
ADC
∆
2
−
1
1
k
ADC
где CADC, CADC < 1 — коэффициент защиты от превышения цифровым сигналом уровня 2
−1;
∆1 — аналоговый эквивалент цифровой «1».
При формировании исходного цифрового радиографического изображения линейным детекто-
ром вектор D будет соответствовать одному детектору линейки и являться основой для построения
строки искомого изображения.
Цифровой сигнал D для устранения влияния искажений, вызванных сцинтиллятором, поступа-
ет на вход процедуры оценки (восстановления) исходного (идеального) цифрового сигнала DI, т.е.
сигнала I(t), подвергнутого интегрированию и дискретизации и оцифровке согласно (6), (7).
Отметим, что с целью возможности сравнения D и DI необходимо использовать один и тот же
аналоговый эквивалент цифровой «1».
2. ВОССТАНОВЛЕНИЕ СИГНАЛА ПОСЛЕ СЦИНТИЛЛЯТОРА
2.1. Формулировка задачи восстановления и ее особенности
Сформулируем задачу восстановления применительно к рассматриваемому случаю. Пусть
имеется сигнал I(t), который изменяется во времени на интервале от 0 до tmax. Сигнал I(t) транс-
формируется (искажается) в сцинтилляционном детекторе функцией h(t), упомянутой выше, в
соответствии с преобразованием (5) и превращается в сигнал J(t).
Уравнение (7) относительно неизвестной функции I(t) по классификации типов интегральных
уравнений относится к линейному интегральному уравнению Вольтерра первого рода типа «сверт-
ки» [18, 19]. Следует отметить, что в случае ассоциации с понятием «исходный сигнал» функции
P(t) (см. формулу (1)), интегральное уравнение (5) превратится в нелинейное интегральное урав-
нение Вольтерра первого рода типа «свертки» [20].
Оценка функции I(t) по измеренной зависимости J(t) и является целью задачи восстановления
в общем случае.
В соответствии с (5), (7):
j
Jd ≈C ∑h(t −s )I(s )∆T, j
=1...
(8)
j
E
j
i
i
jmax
i =1
Существует множество подходов к решению уравнений вида (5), (8), например, [18—21].
Рассмотрим один из них с учетом особенностей формирования и трансформации исходных сиг-
налов.
Для цифровых сигналов система уравнений (8) в матричном виде выглядит следующим
образом:
D = H×DI,
H
=h(
t
−s
),
j=1 ...
j
(9)
ji
j
i
max
Дефектоскопия
№ 6
2022
Восстановление сигналов со сцинтилляционных детекторов
51
Первые две особенности матрицы H связаны с определением передаточной функции h(t) фор-
мулой (2). Первая особенность заключается в том, что h(t) = 0, t < 0, из чего следует, что матрица
H является нижней треугольной матрицей. Вторая особенность обуславливается равенством
max
(0)
(
),
i= j i,
1,...,
j
i
h
=h
t
−s
j=
j
, отсюда вытекает равенствоHjj
f
s
=a
+a
, что существенно
упрощает процесс имитации искажения исходного цифрового сигнала и решения уравнения (9)
относительно DI. Третья особенность вытекает непосредственно из характера функции h(t) и соот-
ношения временных параметров процесса высвечивания в сцинтилляторе и времени интегрирова-
ния ∆t:
h′(t)
<
0;
τ
f
<τ
s
<< t
max
,
∆t << τ
s
,
∆t < τ
f
(10)
Если условие
∆t
<τ
f
не выполняется, то процесс высвечивания можно описать одной экспо-
ненциальной функцией.
Из третьей особенности вытекает существование некоторого момента времени tlim, tlim
<< t
max
,
начиная с которого можно считать, что h(t) = 0. Пусть в момент времени t0 наблюдается идеальный
цифровой сигнал
DI
t
В качестве уравнения нахождения момента времени tlim, отсчитываемого от
0
t0, может быть предложено следующее соотношение:
D(t
)
≈
DI
a
exp
−t
τ
+
a
exp
−t
τ
=ε,
(11)
lim
t
(
f
(
lim
f
)
s
(
lim
s
))
0
где ε — уровень, при котором искажение незначимо. Для цифровых сигналов можно взять, напри-
мер, ε = 1.
Из (11) следует, что значение tlim зависит не только от амплитудных и временных параметров
высвечивания сцинтиллятора τf, af, τs, as, но и амплитуды идеального цифрового сигнала
,
t
DI
0
вызывающего искажение.
Для цифровых сигналов аналогом параметра tlim является параметр jlim. Данный параметр огра-
ничивает суммирование в (8) условием
j-i≤ j
(12)
lim
Из свойств функции h(t) вытекает особенность матрицы H — она является нижнетреугольной:
*
H
0
0
0
0
*
*
H
H
0
0
1
0
*
*
*
H
=
H
H
H
0
,
(13)
2
1
0
*
*
*
*
H
H
H
H
4
3
1
0
*
где
H
=h(k ⋅∆t
),
k=
0, ... ,
j
−1.
k
max
Особенность матрицы (13) сводится к тому, что она полностью определяется одним вектором
*
,
H
что существенно упрощает алгоритм восстановления исходного цифрового сигнала и позво-
ляет реализовать алгоритм обратной свертки.
2.2. Алгоритм восстановления сигнала
Будем считать, что известны параметры, характеризующие сцинтиллятор: τf, af, τs, as. Это озна-
чает, что передаточная функция h(t) полностью определена.
На вход алгоритма поступает дискретизированный цифровой сигнал D, являющийся строкой
цифрового изображения.
Сделаем несколько предварительных замечаний для построения алгоритма восстановления.
Пусть до момента времени t0 светового потока на входе фотоприемника не наблюдалось, то есть
источник излучения был выключен либо полностью экранирован. С моментом времени t0 связыва-
ется появление и регистрация ионизирующего излучения, t0 соотносится с j0 элементом вектора D.
Очевидно, что во все элементы сигнала D, начиная с номера j0 + 1 и заканчивая j0 + jlim, сцинтилля-
тором вносятся значимые вклады от исходного сигнала, т.е. сигнала
DI
Согласно (5) и (8), вос-
j0
становленный сигнал в точке j0 находится по формуле:
Дефектоскопия
№ 6
2022
52
С.П. Осипов, С.А. Щетинкин, Е.Ю. Усачёв и др.
D
∆t
D
∆t
*
j
0
j
0
R
≈
=
j0
*
(14)
a
+a
H
f
s
0
Откорректируем искажение, вносимое сигналом в элементе j0 в последующие элементы:
t
−t
t
−t
*
*
j
j0
j
j0
*
*
D
j
=
D
j
−
R
j
a
f
exp−
+a
s
exp−
∆t
=
D
j
−
R
j
H
j- j
,
j
=
j
0
+1, ... ,
j
0
+
j
lim
(15)
0
0
0
τ
τ
f
s
Корректировка продолжается до элемента с номером j0 + jlim, вклад в который меньше заданно-
го заранее значения ε ( см. (11)).
Если выполняется условие
*
D
j
+1
>ε,
(16)
0
то формула для вычисления значения восстановленного сигнала в элементе с номером j0 + 1 имеет
вид:
*
D
j
+1
∆t
*
0
R
≈
(17)
j
0
+1
*
H
0
По формуле, аналогичной (15), находим искажения, вносимые в элементы частично откоррек-
тированного сигнала, последующие за элементом j0 = 1, вплоть до элемента с номером j0 + 1 + jlim.
Следует отметить, что анализируется уже откорректированный сигнал со своим значением jlim.
*
Если (16) не выполняется, т.е.
1
,
j
D
+
≤ε
то осуществляется поиск дальнейшего по шкале вре-
0
мени ненулевого значения сигнала.
Описанная процедура восстановления продолжается вплоть до достижения последнего эле-
мента цифрового сигнала.
Алгоритм, базирующийся на вышесказанных замечаниях и формулах (14)—(17), представляет
собой обратный ход в методе Гаусса решения системы линейных алгебраических уравнений (9).
Следует отметить, что формально в описанном алгоритме будет фигурировать не вся матрица
H (9), а только ее левый верхний минор A — jlim строк и jlim столбцов. Для такой матрицы A суще-
ствует и матрица A-1. Причем матрица A-1 обладает свойством (13), т.е. определяется первым
-1
столбцом
B
=
(
A
)1
Описанное свойство на практике существенно упрощает процесс восстановления исходного
сигнала, сводя его к использованию обратной свертки (фильтра), где в качестве фильтра F высту-
пает вектор-матрица B с необходимой корректировкой:
B
j
F
=
,
j
=
1
j
j
jlim
lim
(18)
Bi
∑
i =1
Формула для оценки восстановленного сигнала DI* в этом случае имеет вид:
0,
t
<
0∨
t
>
t
j
j
max
*
j
DI
=
F
,
j- i
j
(19)
j
j-i
≤
lim
D
⋅
i
∑
i=0
0,
j- i
>
j
lim
Алгоритм, базирующийся на преобразовании (19), отличается высокой точностью и произ-
водительностью. При малых значениях ∆t остаточное отклонение восстановленного сигнала
от исходного сигнала мало, что следует из положений, использованных при построении алго-
ритма.
Представленный выше алгоритм базируется на информации о параметрах сцинтиллятора τf, τs,
af, as. Имеющиеся в научной литературе данные по этим параметрам могут сильно отличаться друг
от друга. Отсюда вытекает необходимость оценки τf, τs, af, as непосредственно в ходе физического
эксперимента либо в процессе формирования цифровых радиографических изображений. Ниже
обсудим возможные подходы к определению параметров сцинтиллятора τf, τs, af, as.
Дефектоскопия
№ 6
2022
Восстановление сигналов со сцинтилляционных детекторов
53
3. АЛГОРИТМ ОЦЕНКИ ПАРАМЕТРОВ БЫСТРОЙ И МЕДЛЕННОЙ СОСТАВЛЯЮЩИХ
ВЫСВЕЧИВАНИЯ В СЦИНТИЛЛЯЦИОННОМ ДЕТЕКТОРЕ
3.1. Предпосылки разработки алгоритма
Существует физический фактор, позволяющий разработать высокоточный и высокопроизво-
дительный алгоритм оценки временных и амплитудных параметров сцинтилляторов непосред-
ственно в результате физического эксперимента. Упомянутый фактор связан с тем, что временные
параметры быстрой τf и медленной τs составляющих процесса высвечивания в сцинтилляционном
детекторе, как правило, существенно отличаются друг от друга [22, 23], т.е. τf << τs.
Для оценки искомых параметров сцинтиллятора необходимо использовать специальным
образом организованные сигналы, например, ультракороткий или прямоугольный сигналы с
высокой амплитудой на нулевом фоне [4]. Для ультракороткого импульса его длительность timp
должна быть существенно меньше τf, а время измерения ∆t — много меньше timp. Ультракороткий
импульс достаточно сложно физически реализуем. Прямоугольный высокоамплитудный импульс
может быть реализован с помощью сканирования металлической пластины с прорезанной
щелью, причем толщина пластины должна превышать проникающую способность для исследу-
емой максимальной энергии рентгеновского излучения. При этом можно использовать достаточ-
но широкую щель. Предполагается, что информативной частью окажется убывающая часть
цифрового сигнала.
3.2. Аналитическое обоснование алгоритма
Пусть идеальный сигнал I(t) представляет собой прямоугольный высокоамплитудный импульс,
описываемый выражением:
1,
t
≤
t
≤
t
+t
0
0
imp
I t)
=
I
,
(20)
0
0,
(
t
<
t
)
∨
(
t
+t
<
t
<
t
)
0
0
imp
max
где t0, timp — момент прихода и длительность импульса; tmax — длительность сигнала.
После подстановки (20) в (5) и отсчета времени от t0 получим:
t
a
exp
( (
t
−
p
)
τ
)
+
a
exp
( (
t
−
p
)
τ
d
p,
0<
t
≤
t
∫
(
f
f
s
s
))
imp
0
J t)
=
C
I
⋅
(21)
E
0
t
imp
(
a
exp
( (
t
−
p
)
τ
)
+
a
exp
( (
t
−
p
)
τ
))
d
p,
t
<
t
<
t
−t
∫
f
f
s
s
imp
max
0
0
Введем переменную x = t - timp и перепишем вторую часть сигнала Ja(t), т.е. сигнал для
t>t
imp
,
в следующем виде:
C
I
k
⋅
exp
(
−
x
τ
)
+
k
exp
(
−
x
τ
J x)
E
0
(
f
f
s
s
))
a
(22)
W x)
=
=
,
J
(x
=
0)
J
(x
=
0)
a
a
где kf, ks — коэффициенты для быстрой и медленной составляющих высвечивания, не зависящие
от времени t. Формулы для вычисления kf и ks имеют вид:
t
imp
t
imp
k
=a
exp
−x τ
d
p
,
k
=a
exp
(
−x τ
)
d
p
(23)
f
f
∫
(
f
)
s
s
∫
s
0
0
Заметим, что в соответствии с организацией тестового прямоугольного сигнала в ближней
области к импульсу заметно влияние быстрой составляющей, а в дальней — только медленной.
За основу для разработки алгоритма оценки временных и амплитудных параметров сцинтил-
лятора выберем выражение (22) с коэффициентами (23).
Дефектоскопия
№ 6
2022
54
С.П. Осипов, С.А. Щетинкин, Е.Ю. Усачёв и др.
3.3. Этапы алгоритма оценки параметров сцинтиллятора
На вход алгоритма поступает цифровой сигнал D. Известны следующие параметры схемы
контроля: время интегрирования ∆t; скорость сканирования V; ширина щели, формирующей тесто-
вый цифровой сигнал L; уровень цифрового сигнала при измерениях без объекта контроля D0.
Ожидаемая длительность импульса timp оценивается исходя из ширины щели L и скорости скани-
рования V: timp = L/V. Альтернативой формированию цифрового сигнала щелью может быть его
формирование при оптимальном устойчивом режиме работы рентгеновского аппарата и выключе-
нием его (рентгеновского аппарата) через строго заданный промежуток времени timp [4]. В обеих
вариантах формирования тестовых прямоугольных сигналов ожидаемое или заданное значения timp
и D0 корректируются экспериментально.
Обсудим этапы алгоритма.
1. Находим участок изменения координаты j от j1 до j2, на котором значения сигнала D близки
к константе (максимуму) и начиная j = j2 + 1 элементы сигнала убывают. Формула для нахождения
откорректированного значения сигнала D0 имеет вид:
j2
∑
Dj
j= j
1
D
=
(24)
0
j
2
-
j
1
Значение j = j2 соответствует
t
=t
j
2
imp
2. Определяем значение j, j = j3 по экспериментальному представлению формулы (11):
D
=
min
{
j,
j
< j≤ j
D
=ε
}
(25)
j
3
2
max
j
3. Формируем экспериментальную дискретизацию функции затухания W0, W1,..,
W
с помо-
i
lim
щью выражения:
D
i+ j2
W
=
,
i
=
0, ... ,i
,
i
=
j
−
j
(26)
i
lim lim
3
2
D
0
4. Значения элементов вектора W прологарифмируем:
LW
=
lnW
,
i=
0, ... ,i
(27)
i
i
lim
5. Находим координату iL, 0 < i < ilim, для которой линия, соединяющая точки
i
,lnD
и
(
L
i
L
)
(
i
,lnD
)
,
близка к прямой. Координата iL находится движением от точки ilim в сторону умень-
lim
i
lim
шения с
подсчетом текущих параметров линии методом наименьших квадратов. Если сигнал в
*
точке i3 достаточно близок к нулю, то координата i уменьшается до некоторого значения
i=i
,
3
начиная с которого цифровой сигнал начинает возрастать при уменьшении i. После уменьшения
значение координаты i на несколько точек вычисляются текущие параметры линии и величина
невязки Δ, отнесенная к общему числу точек n, вовлеченных в процедуру обработки. Процедура
*
вычислений прерывается в точке
i=i
, если отношение ∆/n, начиная с этой точки, начинает резко
2
и монотонно возрастать.
*
*
6. В соответствии с подходом, отмеченным выше, на интервале
j
,
j
находятся параметры
2
3
медленной составляющей — as и τs.
Для этого реализуем несколько шагов.
Введем функцию:
Flin(x,a,b) = -ax +b
(28)
Зададим минимизируемую функцию:
*
2
i3
*
Fs a,b)
=
Flin
i- i
∆t a,b
−
LW
(29)
∑
(
((
2
)
)
i
)
*
i=i2
Введем начальные значения a и b, а затем найдем минимальное значение функции Fs(a,b) по
параметрам a и b. Функция Fs(a,b) достигает минимум для значений параметров a = as, b = bs ,
если выполняется равенство:
Дефектоскопия
№ 6
2022
Восстановление сигналов со сцинтилляционных детекторов
55
min
a,b
Fs a,b)
= Fs(as,bs).
(30)
Искомые параметры τs, as находятся по формулам:
1
exp(bs)
τ
=
,
a
=
(31)
s
s
t
imp
as
exp
(
−
x
τ
)
d
p
∫
s
0
7. Оценивается вклад быстрой составляющей в цифровой сигнал:
(
i- i
)
∆t
2
*
Df
=
D
−a
×exp
−
,i
=
i
,...,i
(32)
i
i
s
2
2
τ
s
8. По аналогии с пунктом 3 формируем экспериментальную дискретизацию функции затуха-
ния для быстрой составляющей Wf0, Wf1, ... , Wf с помощью выражения:*
i
2
*
Wf
= Df
Df
,
i=
0, ... ,i
(33)
i
i
0
2
9. Значения элементов вектора Wf прологарифмируем:
*
LWf
=
lnWf
,
i=
0, ... ,
i
(34)
i
i
2
**
**
*
10. Для нахождения соответствующего участка
i
, ... ,i
,i
<i
и для оценки параметров af и τf
2
2
2
2
воспользуемся подходом, описанным в пунктах 3—5, с той лишь разницей, что координата i не
уменьшается, а увеличивается, и анализируется сигнал Df.
**
11. На интервале
j
2
,
j
2
находятся параметры медленной составляющей — as и τs.
Для этого воспользуемся введенной ранее функцией (28).
Аналогично (29) зададим минимизируемую функцию:
**
i
2
2
*
(35)
Ff a,b)
=
∑
(
Flin
((
i- i
2
)
∆
t a,
b
)
−
LWf
i
)
i=0
Введем начальные значения a и b, а затем найдем минимальное значение функции Ff (a,b) по
параметрам a и b. Функция Ff (a,b) минимальна для значений параметров a = af , b = bf , если
выполняется равенство:
min
Ff a,b)
= Ff
(
af
,
bf
).
(36)
a
,b
Параметры быстрой составляющей высвечивания τf, af находятся по формуле, аналогичной
(31):
t
imp
1
τ
=
,
a
=
exp(bf
)
exp
(
−
x
τ
)
d
p
(37)
f
f
f
∫
af
0
Предложенный алгоритм реализован в системе для математических вычислений MathCad.
Отметим, что поиск линейных участков может быть осуществлен и на основе визуального
анализа соответствующих участков графиков.
4. РЕЗУЛЬТАТЫ ОБРАБОТКИ ДАННЫХ ВЫЧИСЛИТЕЛЬНЫХ ЭКСПЕРИМЕНТОВ
Для апробации предложенных выше алгоритмов восстановления сигнала после сцинтилля-
ционного детектора и оценки амплитудных и временных параметров быстрой и медленной
составляющих затухания сигнала в сцинтилляционном детекторе был проведен ряд вычисли-
тельных экспериментов. В качестве сцинтиллятора был выбран GSO. Для данного материала в
[23] приведены следующие значения временных параметров: τf = 65,2 нс; τs = 605,8 нс и соот-
ветствующий вклад в интегральную интенсивность If = 86 %, Is = 14 %. В соответствии с (5),
значения амплитудных параметров af и as равны 6,104×106 и 9,937×105 соответственно.
Дефектоскопия
№ 6
2022
56
С.П. Осипов, С.А. Щетинкин, Е.Ю. Усачёв и др.
Сигнал со сцинтилляционного детектора оцифровывался АЦП с разрядностью kADC = 16.
Время регистрации ∆t = 25 нс. Коэффициент защиты цифрового сигнала от превышения макси-
мального значения CADC = 0,8. Для этих данных D0 = 52428.
В табл. 1 приведены результаты расчетов по формуле (11) отношения tlim/∆t в целых числах для
приведенных выше условий и диапазона изменения ε от 1 до 10. Отмеченное отношение является
количеством элементов jlim фильтра F в обратной свертке (19) для D0 = 52428. Для обоснования
необходимости использования фильтра с высоким количеством элементов в обратной свертке в
табл. 1 также указаны отношения tlim/τs.
Таблица
1
Результаты расчетов отношения tlim/∆t (11) для D0 = 52428
ε
1
2
3
4
5
6
7
8
9
10
tlim/∆t
269
246
240
233
223
218
215
210
208
205
tlim/τs
8,87
8,12
7,93
7,69
7,37
7,2
7,1
6,93
6,88
6,78
Приведенные в табл. 1 данные свидетельствуют о необходимости учета медленной компонен-
ты процесса высвечивания энергии фотонов, поглощенных в сцинтилляционном детекторе.
4.1. Результаты восстановления сигнала после сцинтилляционного детектора
Выбор исходного сигнала для тестирования алгоритма восстановления сигнала после сцин-
тиллятора зависит от задач, стоящих перед системами ЦР и РКТ. В качестве обязательных тре-
бований к такого рода сигналам можно придерживаться требования перекрытия возможного
диапазона изменения сигнала и необходимости обнаружения локальных сигналов малого кон-
траста [16, 17, 24].
В качестве объектов для симулирования искажений, вносимых сцинтиллятором в исходный
сигнал, были выбраны ступенчатые и клиновидные объекты с локальным включениями в форме
кубовидных пор. При моделировании исходных объектов придерживались подхода, изложенного
в работах [16, 17].
Отметим, что информативному сигналу предшествовало полное экранирование и часть сигна-
ла, необходимая для калибровки по белому.
На рис. 1 приведены идеальные сигналы, их искажения сцинтилляционным детектором и
результаты реконструкции для клиновидных и ступенчатых сигналов. Тестовые объекты содержа-
ли также кубические поры размером около 5 % от толщины объекта в месте локализации.
Из сравнения идеальных сигналов на рис. 1а с сигналами, искаженными сцинтилляционным
детектором на рис. 1б, можно сделать вывод о значимом искажении сигналов детектором, особен-
но для участков с резким перепадом сигналов. Наличие малоконтрастных локальных неоднород-
ностей по искаженным сигналам подтверждается, но сделать корректное заключение о форме этих
неоднородностей затруднительно. Восстановленные сигналы близки к исходным, незначительные
остаточные искажения при реконструкции наблюдаются в местах с резким и значительным изме-
нением амплитуды сигналов.
4.2. Результаты оценки параметров быстрой и медленной составляющих высвечивания
сцинтиллятора
В алгоритме оценки параметров быстрой и медленной составляющих высвечивания сцинтил-
ляционного детектора отмечена необходимость исследования затухающей части высокоамплитуд-
ных прямоугольных сигналов. Анализ данных, приведенных на рис. 1, подтверждает существова-
ние выраженной ниспадающей части сигнала.
На рис. 2 приведено изображение прямоугольного сигнала, искаженного сцинтилляционным
детектором. При симуляции использовались параметры из предыдущего пункта, исключение
∆t = 10 нм. Сигнал зашумлен мультипликативным шумом (распределение Гаусса, σ = 0,01).
Дефектоскопия
№ 6
2022
Восстановление сигналов со сцинтилляционных детекторов
57
Клиновидный объект
Ступенчатый объект
DI
а
DI
60 000
60 000
50 000
50 000
40 000
40 000
30 000
30 000
20 000
20 000
10 000
10 000
0
0
0
800
1600
2400
3200
t/∆t
0
800
1600
2400
3200
t/∆t
б
D
D
60 000
60 000
50 000
50 000
40 000
40 000
30 000
30 000
20 000
20 000
10 000
10 000
0
0
0
800
1600
2400
3200
t/∆t
0
800
1600
2400
3200
t/∆t
в
R
R
60 000
60 000
50 000
50 000
40 000
40 000
30 000
30 000
20 000
20 000
10 000
10 000
0
0
0
800
1600
2400
3200
t/∆t
0
800
1600
2400
3200
t/∆t
Рис. 1. Сравнение сигналов для клиновидного и ступенчатого объектов:
а — идеальные сигналы; б — сигналы, искаженные сцинтилляционным детектором; в — восстановленные сигналы.
D
D
а
б
60 000
60 000
50 000
50 000
40 000
40 000
30 000
30 000
20 000
20 000
10 000
10 000
0
0
0
800
1600
2400
3200
t/∆t
3000
3100
3200
t/∆t
Рис. 2. Изображение прямоугольного сигнала, искаженного сцинтилляционным детектором, (a) и трансформация иде-
ального сигнала (б).
Дефектоскопия
№ 6
2022
58
С.П. Осипов, С.А. Щетинкин, Е.Ю. Усачёв и др.
С целью иллюстрации возможности применения исследуемого алгоритма на рис. 2 продемонстри-
рована затухающая часть сигнала.
Из визуального анализа ниспадающей части искаженного сцинтиллятором сигнала можно сде-
лать вывод о возможности использования этой информации для оценки параметров быстрой и
медленной составляющих высвечивания.
Для проверки качества алгоритма оценки временных параметров был проведена серия вычис-
лительных экспериментов по моделированию процесса трансформации идеального сигнала сцин-
тилляционным детектором с последующей оценкой временных параметров для быстрой τf и мед-
ленной τs составляющих процесса высвечивания сцинтилляций. При этом варьировались уровни
шумов сигналов D(t).
В табл. 2 приведена зависимость оценок параметров τf и τs от уровней относительных
шумов δ.
Таблица
2
Оценки параметров τf и τs от уровней относительных шумов δ
δ, %
Параметры
Исходные
0
0,2
0,4
0,6
0,8
1
2
af
6,104×106
6,104×106
6,094×106
6,104×106
6,105×106
6,098×106
6,078×106
6,134×106
τf, нс
65,2
65,2
64,8
64,6
64,8
64,8
64,8
64,7
as
9,937×105
9,937×105
9,957×105
9,972×105
9,958×105
9,951×105
9,947×105
9,990×105
τs, нс
605,8
605,8
605,4
605,3
605,4
605,4
605,4
605,3
Из анализа данных, представленных в табл. 1, можно сделать вывод о достаточно высокой
устойчивости предлагаемого алгоритма оценки параметров τf и τs к шумам. Оценки погрешности
исследуемых параметров близки к экспериментальным погрешностям измерения τf и τs, указанным
в [23]. Вероятнее всего высокий уровень оценки погрешности τf по сравнению с погрешностью τs
обусловлен для рассматриваемого примера ограниченностью интервала времени, по которому и
оценивался параметр быстрой составляющей. Очевидно, что с увеличением времени измерения
точность оценки параметров быстрой составляющей высвечивания и чувствительность к шумам
будет резко возрастать.
5. ОБСУЖДЕНИЕ И РЕКОМЕНДАЦИИ
5.1. Общие замечания
Для нахождения параметров Bf и τf может быть использована начальная часть цифрового сиг-
нала, соответствующая анализируемому импульсу. В этом случае медленная составляющая долж-
на быть откорректирована.
Отметим, что методы многопараметрической минимизации (четыре параметра) недостаточно
эффективны применительно к рассматриваемой задаче по сравнению с предлагаемым алгоритмом,
который отличается простой, ясностью и легкостью программной реализации.
Разумеется, что качество алгоритмов, предложенных в разделах 2 и 3, определяется помимо
параметров сцинтиллятора также формой и длительностью сигналов, амплитудами и уровнями
шумов цифровых сигналов, вдобавок временем измерения (интегрирования) ∆t.
Заметим, что коррекция послесвечения применительно к системам ЦР и РКТ должна быть
осуществлена практически в любых случаях, за исключением использования сверхбыстрых
сцинтилляторов и контроля объектов без резких перепадов по толщине в длинах свободного
пробега.
Предложенные выше алгоритмы реализованы в системе для математических вычислений
MathCad. Данный выбор обусловлен максимальной близостью языков математики и MathCad
[16, 17]. Программы, написанные в системе MathCad, отличаются понятностью, адаптивностью и
простотой отладки.
Дефектоскопия
№ 6
2022
Восстановление сигналов со сцинтилляционных детекторов
59
Разработанные алгоритмы и программы их реализующие естественным образом объединяются
с существующими математическими и имитационными моделями систем цифровой радиографии
[16, 17] и рентгеновской компьютерной томографии [25].
Существует ряд вопросов, связанных с рассматриваемыми задачами, для решения которых
можно применить разработанные алгоритмы. Остановимся подробнее на одном из них.
5.2. Устойчивость алгоритма оценки параметров высвечивания
Выше изложен алгоритм оценки амплитудных и временных параметров быстрой и медленной
составляющих высвечивания af, as, τf, τs. Одной из главных характеристик любой задачи измерений
параметров является оценка устойчивости соответствующего алгоритма обработки информации.
Определим понятие устойчивости алгоритма применительно к рассматриваемой задаче. Пусть
имеется вектор оцениваемых параметров
P=
(
a
f
τ
f
a
s
τ
s
).
В процессе измерений формиру-
ется информационный сигнал D, который также является вектором, но с существенно большей
размерностью. Между векторами P и D существует связь, описываемая некоторым оператором F:
D = F(P).
(30)
Алгоритм оценки амплитудных и временных параметров высвечивания реализует оператор
-1
-
F
, который близок к обратному оператору
1.
F
В результате реализации алгоритма применитель-
но к экспериментальным данным получим преобразование:
−1
−1
*
P
=
F
(
D
)
=
F
(
F P)
)
(31)
*
В формуле (31) вектор
P
— это вектор оценок исследуемых параметров высвечивания
*
*
*
*
a
τ
a
τ
(
à
f
s
s
)
Понятие устойчивости алгоритма связывают с его реакцией на малые возмущения исходных
данных [21]. Применительно к рассматриваемой задаче можно принять, что алгоритм определения
параметра p называется устойчивым, если малое приращение значения исходного параметра ∆p
*
приводит к малому приращению оценки исследуемого параметра
∆p
и выполняется условие:
*
∆p
∆p ≈ 1.
(32)
На практике предпочтительным является проведение исследований при вариации параметров
в относительных величинах, поэтому (32) трансформируется в ограничение:
*
δp
δp ≈ 1.
(33)
Проведение исследование устойчивости алгоритма для всей совокупности сцинтилляторов в
рамках одной статьи не представляется возможным и в этом нет особой необходимости.
Проанализируем устойчивость для приведенного в пункте 4.2 примера.
Результаты моделирования показали близость рассматриваемого отношения для всех информа-
тивных параметров к единице при изменении параметров до 10 и более процентов. Это связано с
тем, что для уровней шумов в исходных сигналах близких к нулю автоподстройка алгоритма
устремляет систематическую погрешность оценки исследуемых параметров к нулю. Расхождение
же связано исключительно с цифровым эквивалентом единицы и смещением реального и ожидае-
мого времени положения максимума для прямоугольного сигнала.
ЗАКЛЮЧЕНИЕ
Предложена математическая модель искажения исходного радиометрического сигнала сцин-
тилляционным детектором для случая регистрации непрерывного рентгеновского излучения.
Модель учитывает распределение интенсивности излучения по аналоговому сигналу, время
регистрации излучения, временные и амплитудные параметры быстрой и медленной составля-
ющих процесса высвечивания в сцинтилляционном детекторе. Разработан способ восстановле-
ния сигнала после сцинтиллятора, основанный на предварительной информации о параметрах
сцинтилляционного детектора. Предложен алгоритм оценки амплитудных и временных параме-
тров быстрой и медленной составляющих процесса высвечивания в сцинтилляционном детекто-
Дефектоскопия
№ 6
2022
60
С.П. Осипов, С.А. Щетинкин, Е.Ю. Усачёв и др.
ре. Проведена серия вычислительных экспериментов и на основе анализа ее результатов проде-
монстрирована эффективность разработанных алгоритмов. Разработанные алгоритмы восста-
новления сигналов, искаженных сцинтилляционным детектором, и оценки параметров быстрой
и медленной составляющих процесса высвечивания могут быть усовершенствованы в части
детального описания начального этапа высвечивания и представления затухания суммой из трех
и более экспоненциальных функций.
СПИСОК ЛИТЕРАТУРЫ
1. Kumar V., Luo Z. A review on x-ray excited emission decay dynamics in inorganic scintillator materials
2. Lu L., Sun M., Wu T., Lu Q., Chen B., Huang B. All-inorganic perovskite nanocrystals: new generation
scintillators for high-resolution X-ray imaging // Nanoscale Advances. 2022. V. 4. P. 680—696. https://doi.
org/10.1039/D1NA00815C
3. Yasar F., Kilin M., Dehdashti S., Yu Z., Ma Z., Wang Z. Spatially resolved x-ray detection with photonic
crystal scintillators // Journal of Applied Physics. 2021. V. 130. No. 4. No. article 043101. https://doi.
org/10.1063/5.0050380
4. Nagarkar V.V., Vasile S., Gothoskar P., Gordon J.S., Gupta T.K. CCD based high resolution non-
destructive testing system for industrial applications // Applied Radiation and Isotopes. 1997. V. 48.
5. Panetta D. Advances in X-ray detectors for clinical and preclinical Computed Tomography // Nuclear
Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and
6. Eilbert R.F. X-ray technologies // Aspects of explosives detection. Amsterdam: Elsevier,
2009.
7. Foster C., Wu Y., Stand L., Koschan M., Melcher C.L. Effect of lithium codopant concentration on the
luminescence properties of (Lu0. 75Y0. 25) 3Al5O12: Pr3+ single crystals: Before and after air annealing //
8. Marshall M.S., Kenesei P., Marton Z., Sosa C., Brecher C., Wart M., Miller S., Singh B., Miceli A.,
Nagarkar V.V. Advances in high-resolution ultrafast LuI 3: Ce scintillators for fast timing applications //
TNS.2020.2985213
9 Postupaev V.V. High dynamic range measurement of CdWO4 afterglow // Nuclear Instruments and
Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment.
10. Van der Sar S.J., Brunner S.E., Schaart D.R. Silicon photomultiplier-based scintillation detectors for
photon-counting CT: A feasibility study // Medical Physics. 2021. V. 48. No. 10. P. 6324—6338. https://doi.
org/10.1002/mp.14886
11. Wang X., Dai Y., Zhang Z., Su L., Kou H., Wang Y., Wu A. Optical and scintillation properties of Ce:
Y3Al5O12 single crystal fibers grown by laser heated pedestal growth method // Journal of Rare Earths. 2021.
12. Chen S., Jiang B., Wang Y., Zhu Q., Yang Q., Ma W., Zang G., Zhang L. Fabrication of Ce-doped
(Gd2Y) Al5O12/Y3Al5O12 composite-phase scintillation ceramic // Journal of Rare Earths. 2019. V. 37. No. 9.
13. Nagorny S. Novel Cs2HfCl6 Crystal scintillator: recent progress and perspectives // Physics. 2021.
14. Wu T., Wang L., Shi Y., Xu T., Wang H., Fang J., Ni. J., He H., Wang Ch., Wan B., Ding D., Zhou Zh.
Liu Q., Li Q., Ju J., Huang X., Shichalin O., Papynov E.K. Fast (Ce, Gd)3Ga2Al3O12 scintillators grown by the
optical floating zone method // Crystal Growth & Design. 2021. V. 22. No. 1. P. 180—190. https://doi.
org/10.1021/acs.cgd.1c00779
15. Behling R. X-ray sources: 125 years of developments of this intriguing technology // Physica Medica.
16. 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.
лов С.В., Осипов О.С., Усачёв Е.Ю. Имитационное моделирование в цифровой радиографии с учетом
пространственных форм объектов контроля // Дефектоскопия. 2020. No. 8. С. 35—48.]
17. 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. No. 2. С. 43—55.]
Дефектоскопия
№ 6
2022
Восстановление сигналов со сцинтилляционных детекторов
61
18. Jaabar S.M., Hussain A.H. A Move recent review of the integral equations and their applications //
Journal of Physics: Conference Series. — IOP Publishing, 2021. V. 1818. No. 1. No. article 012170. https://
doi.org/10.1088/1742-6596/1818/1/012170
19. Karchevsky A.L. Solution of the convolution type Volterra integral equations of the first kind by the
quadrature-sum method // Journal of Applied and Industrial Mathematics. 2020. V. 14. No. 3. P. 503—512.
20. Alvarez E., Lizama C. Attractivity for functional Volterra integral equations of convolution type //
cam.2016.01.048
21. Lubich C. Convolution quadrature revisited // BIT Numerical Mathematics. 2004. V. 44. No. 3.
22. Usui Y., Nakauchi D., Kawano N., Okada G., Kawaguchi N., Yanagida T. Scintillation and optical
properties of Sn-doped Ga2O3 single crystals // Journal of Physics and Chemistry of Solids. 2018. V. 117.
23. Nassalski A., Kapusta M., Batsch T., Wolski D., Mockel D., Enghardt W., Moszynski M. Comparative
study of scintillators for PET/CT detectors / IEEE Nuclear Science Symposium Conference Record. 2005 //
24. American national standard for determination of the imaging performance of X-ray and gamma-ray
systems for cargo and vehicle security screening. Standard ANSI N42.46, 2008.
25. Osipov S.P., Yadrenkin I.G., Chakhlov S.V., Osipov O.S., Usachev E.Yu., Manushkin A.A. Calculation
model of X-ray computed tomography with density assessment function // Russian Journal of Nondestructive
Ядренкин И.Г., Чахлов С.В., Осипов О.С., Усачёв Е.Ю., Манушкин А.А. Вычислительная модель рентге-
новской компьютерной томографии с функцией оценки плотности // Дефектоскопия. 2021. № 3.
С. 37—52.]
Дефектоскопия
№ 6
2022
