Сцинтилляционная гамма спектрометрия

 Рассмотрим основные проблемы, которые приходится решать при обработке сцинтилляционных гамма спектров и методы их решения, реализованные в программном пакете «СПЕКТР»:

Энергетическое разрешение

 
 
Рис. 1. Спектр Ra-226 с продуктами распада.
 
Во времена царя Гороха, когда считать приходилось пользуясь логарифмическими линейками и арифмометрами, активность нуклидов в пробах измеренных на ППД спектрометрах определялась по отношению площадей отдельных пиков в спектре пробы и спектре эталона. В сцинтилляционной спектрометрии выделить отдельный пик вследствие плохого энергетического разрешения спектрометра можно только в крайне редких искусственных случаях. На рис. 1 представлен спектр Ra-226 с дочерними продуктами распада в относительно равновесном состоянии измеренный на детекторе из особо чистого германия (нижний спектр) и спектр той же пробы измеренный на сцинтилляционном детекторе NaI(Tl) (верхний спектр). Разрешение ППД спектрометра - 2 кэВ по линии 1332 кэВ, разрешение сцинтиллятора - 8% по линии 661 кэВ – средние значения для большинства аналогичных детекторов.
 
Спектры проб содержащих естественно радиоактивные нуклиды, кроме Ra-226 содержат ещё не менее богатые пиками цепочки распада Th-232 (11 нуклидов в цепочке распада и более 150 значимых гамма линий, не считая линий сложения каскадных и линий утечки), цепочку U-235 (13 нуклидов в цепочке распада и около 200 заметных гамма линий) и цепочку U-238 – Th-230 (6 нуклидов и масса мелких и не очень гамма линий – более 250).
 
Невозможность для реальных сцинтилляционных спектров выделения каких-либо пиков (не говоря уже о «чистых» пиках отдельных нуклидов) заставляла решать систему уравнений сформированных из интегральных энергетических окон. Невысокая скорость расчётов на логарифмической линейке ограничивала количество окон и не позволяла оценивать какие-либо дополнительные параметры решаемой задачи, кроме активности нуклидов.
 

Энергетический дрейф

 
И если в ППД спектрометрии процедура оценки активности нуклида по площади пика устойчива к энергетическому дрейфу спектрометра (мы можем ошибиться с идентификацией пика, но площадь его останется неизменной при любом дрейфе спектрометра), оконный метод в сцинтилляционной спектрометрии весьма чувствителен к энергетическому дрейфу спектрометра. А сцинтилляционные тракты, к сожалению, дрейфуют – то есть коэффициент преобразования сцинтилляционных детекторов не стабилен, а изменяется с течением времени. Дрейф зависит от времени работы тракта, температуры, интенсивности регистрируемого излучения и стабильности питания.
 
Всевозможные схемы стабилизации работают плохо. Почему? Да потому что невозможно без качественной полной обработки спектра решить задачу стабилизации. Сплощь и рядом в схемах стабилизации используется обратная связь (коэффициент усиления упал – поднимем коэффициент усиления и наоборот). Но любая обратная связь с инерционными элементами (а безинерционных элементов в природе не существует) приводит к осцилляциям (те, кто помнят о преобразовании Лапласа и передаточной функции, спорить не будут), а значит, разрешение портится. А нет никакой нужды использовать обратную связь, когда у нас под рукой компьютер – можно просто пересчитывать спектр, либо добавлять при регистрации импульса не единичку в один канал, а, например, 0.8 в один канал, и 0.2 в соседний, в соответствии со значениями параметров дрейфа.

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

При попытке использовать для коррекции только один пик, три параметра мы никак определитиь не можем. Мы и один то параметр определяем очень плохо. Почему? Предположим, мы используем светодиод генеририующий ипмульсы определённой амплитуды закреплённый у ФЭУ. В этом случае мы можем, конечно, генерировать импульсы любой интенсивности и проблема малой статистики снимается, но остаются проблемы стабильности импульсов от светодиода, проблема коррекции мёртвого времени и проблема только одного реперного пика. Если мы используем реперный источник – он портит спектр, да и наложение измеряемого спектра на репер искажает стабилизацию. Предположим мы используем пики самой пробы или фоновый калий-40. Но при этом различный нуклидный состав пробы может как угодно смещать стабилизацию, поскольку излучение других нуклидов (в том числе и K-40 – его положение в спектре отличается от положения фонового K-40) накладывается на фоновый K-40. Не поможет и стабилизация по пику 2614.5 кэВ Tl-208, поскольку пики сложения каскадных линий от нуклидов цепочки распада Ra-226 накладываются и на этот пик.

Замечательно работает только одна схема стабилизации – полная обработка всего спектра, включающая в качестве определяемых параметров и параметры энергетического дрейфа. В математике, как это может быть и не прискорбно, 2х2=4. Поэтому эффективных оценок (в том числе и оценок параметров дрейфа) можно добиться, только используя всю корректную информацию об объекте корректно.
 
Посмотрим на результаты обработки «сдрейфовавшего» спектра:
  
Рис. 2. Обработка спектра оконным методом без коррекции дрейфа.
 
 
Рис. 3. Обработка спектра стандартным методом с коррекцией дрейфа.
  
 
Наименование объекта исследования: Эталон Th-232. Акт = 977 Бк/кг.
Дата отбора: 01.08.2006 13:48:59
Масса счётного образца: 1638 г
Дата измерения: 10.01.2008 16:25:16
Время измерения: 1800 с
 
Радионуклид Результат исследований (Бк/кг)
Оконный метод.
Результат исследований (Бк/кг)
Стандартный метод.
137Cs 110 ± 210 2 ± 6.3
40K 300 ± 1500 0 ± 59
226Ra 190 ± 340 0 ± 11.6
232Th 620 ± 380 976 ± 76
Невязка: 27.1   (Норма - не более 2)
Энергетич. диапазон: 200 - 2852.24 кэВ.
Дрейф: Лин. 0% (Ку=0.956, Нш=3 кан.),
Нелин. 0%, Разрешения 0%, Zeff +0.0
Невязка: 1.6    (Норма - не более 2)
Энергетич. диапазон: 200 - 2989.75 кэВ.
Дрейф: Лин. 5.1% (Ку=1.005, Нш=1 кан.),
Нелин. 0%, Разрешения -1.67%, Zeff +0.0

 Табл. 1. Результаты обработки

 
На рис. 2 представлена итоговая картинка разложения спектра пробы (в качестве пробы эталон Th-232) оконным методом (7 интегральных окон), на рис. 3 – картинка разложения спектра стандартным для пакета «Спектр» методом и в таблице 1 – результаты обработки тем и другим методом.
 
Из сравнения ясно, что даже небольшой (вполне обычный для сцинтилляционных спектрометров) 5-процентный энергетический дрейф приводит при решении системы методом без коррекции энергетического дрейфа к катастрофическим результатам. Дрейф спектрометра может быть связан с активностью пробы и, в этом случае, не компенсируется предварительной процедурой энергетической калибровки.
 
Процедура решения линейной системы в пакете «Спектр» корректно оценивает погрешности решения и поэтому, хотя результаты решения задачи «оконным методом» никуда не годятся (слишком большие неопределённости оценок активностей), математически они вполне корректны (истинные значения лежат в пределах доверительного интервала) и поэтому не приведут к неверному заключению. Не все программы настолько корректны в оценке погрешностей и в данном случае вполне могли бы выдать, например, результат по Cs-137 110 ± 20 – а такой результат ошибочен и может привести к неверным выводам.
 

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

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

 
Модель спектра при этом выглядит следующим образом:
 
 Здесь:
– спектр фона.
– спектры компонент (нормированные на один Беккерель спектры излучения отдельных нуклидов или цепочек нуклидов, находящихся в состоянии радиоактивного равновесия).
 – амплитуды соответствующих компонент.
и  – параметры энергетического дрейфа спектрометра.
 Множитель .
 
 Z- отклонение условного Z эффективное пробы от Z эффективное эталона использованного при расчёте эффективности регистрации.  
 E(i)- энергия канала i.
- массовый коэффициент поглощения.
- плотность счётного образца. 
- толщина счётного образца (задаётся при установке геометрии).
 Pn(i)- полином (полином используется только в случае предположения о неполном нуклидном составе модели).
 i- канал спектра.
 
 
Свёртка с функцией рассеяния  моделирует ухудшение разрешения спектрометра. В качестве функции рассеяния, за исключением альфа спектрометрии, используется гауссиан . Для альфа спектрометрии используется свёртка с откликом спектрометра на квант излучения - моделью пика .
 
В общем случае предполагаются независимые значения параметров дрейфа для всех компонент. Если все компоненты разложения находятся в единой энергетической шкале, оператор может указать в табличке параметров декомпозиции, что спектры компонент находятся в единой шкале, то есть , ,.
Оператор может исключить из списка параметров  и (или)    и (или)  и , удалив соответствующие галочки в таблице параметров декомпозиции.
 
Решение
Используем метод Гаусса – Ньютона. Начальное значение вектора параметров находим, используя глобальное варьирование выбранных параметров в заданном диапазоне с заданным шагом. Быстродействие современных компьютеров, пока что не позволяет варьировать все параметры с достаточной плотностью сетки, поэтому оператор может задать в параметрах декомпозиции варьирование четырёх основных параметров: Z эффективное пробы, коэффициента усиления, нуля шкалы и разрешения спектрометра. Найденное приближение уточняем, используя ньютоновский процесс.
Здесь:
    – вектор параметров на k –ой итерации;
  – вектор параметров на k+1 –ой итерации.
     – оператор проектирования получаемых значений параметров на ограничения.
    – шаг на k – ой итерации.  
  - матрица Гаусса - Ньютона;
 - матрица весов.   - ковариационная матрица вектора .
 - коэффициент (задается в параметрах декомпозиции) и демпфирующая матрица k -ого шага.
 - ковариационная матрица априорных данных.
 
 
 - априорные данные об активностях компонент.
 - оценка вклада компонент на k -ом шаге,
 - оценка коэффициента усиления тракта.
 - границы МНК оценивания.
 
 
 
Значения активностей на момент начала измерений
 
El, Er - левая и правая энергетические границы декомпозиции, e - энергия;
 - спектр j-ой компоненты (нормированный на один Беккерель спектр излучения нуклида или цепочки нуклидов, находящихся в состоянии радиоактивного равновесия);
 - значения параметров энергетического дрейфа на последней итерации;
 - начальные значения параметров энергетического дрейфа: 
;
- последняя итерация;
 - время измерения;
 - период полураспада нуклида (период полураспада материнского нуклида, обеспечивающего радиоактивное равновесие цепочки нуклидов, в случае цепочки нуклидов).
 
Оценка погрешностей
Мы предполагаем, что спектр   имеет диагональную ковариационную матрицу  ( – единичная матрица,  – вектор дисперсий) причем  для и  для 
 
Здесь:
– количество отсчётов в i – том канале спектра,
 – коэффициент линейной составляющей статистики . То есть спектр имеет составную: пуассоновскую плюс линейную статистику. Для спектра пробы   определяется дифференциальной нелинейностью спектрометра. Для спектра фона   является коэффициентом изменчивости фона, рассчитываемой процедурой статистической обработки результатов измерений фона.
 
Значения в каналах спектра фона полагаем коррелированными (если фон увеличился, то он увеличился, вероятно, во всех каналах).  для . Такая коррелированность сохраняет свойство составной пуассоновской статистики для дисперсии .
 
Здесь:
 
– количество отсчётов в i – том канале спектра фона.
 – коэффициент изменчивости фона. Значение  задаётся в параметрах настройки спектрометра.
В итоге:
 
Здесь
- оценки ковариационной матрицы определяемых параметров на последнем и предпоследнем шаге процедуры решения.
 
и
  - оценки параметров на последнем и предпоследнем шаге процедуры решения.
То есть мы используем оценки погрешностей линеаризованной модели дополненной погрешностями, связанными с обрывом итерационного процесса (погрешности дрейфа спектрометра).
- матрица Гаусса - Ньютона рассчитанная на последнем шаге процедуры решения. Наша процедура решения формирует ковариационную матрицу параметров в матрице обратной матрице Гаусса - Ньютона.
 - матрица весов.  - ковариационная матрица вектора . Считаем её диагональной  (  – единичная матрица,  – вектор дисперсий) причем , где  отношение времени измерения фона и пробы.
 - ковариационная матрица априорных данных.
 
Дисперсии оценок активностей .
 - парциальные (рассчитанные в пределах объединения значимых областей определения коррелированных с i-м компонентом элементарных спектров) невязки для данного нуклида. Значимая область определения - область, где значения элементарного спектра отличны от нуля.
 - диагональный элемент ковариационной матрицы параметров для параметра 
;
- последняя итерация;
 - время измерения;
 - период полураспада нуклида (период полураспада материнского нуклида, обеспечивающего радиоактивное равновесие цепочки нуклидов, в случае цепочки нуклидов).
 - коэффициент линейного энергетического дрейфа ( ) j-ой компоненты рассчитанный как коэффициент линейной регрессии оценённого нелинейного энергетического дрейфа.
 
Итоговая полная погрешность i - ой активности:
       (1)    или
.         (2)
 - дополнительная относительная погрешность - задается как параметр декомпозиции.   связана с погрешностями определения массы пробы, нарушений геометрии и т.д. и т.п.
 - относительная погрешность, рассчитанных при калибровке спектрометра, значений интенсивности компонент разложения и связана с погрешностями аттестации эталонов, погрешностями расчёта элементарных спектров.
 
Поскольку с чиновниками от науки не соскучишься, в зависимости от версии пакета, случайные и систематические погрешности складываются либо линейно (1) в соответствии с рекомендациями Государственного комитета РФ по стандартизации, метрологии и сертификации (МИ 2453-98), либо квадратично (2) в соответствии с рекомендациями Государственного комитета РФ по стандартизации, метрологии и сертификации (МИ 2453-2000). В версиях программы 17.1 и выше складываются квадратично в соответствии с последними рекомендациями.
Оценка погрешности соответствует 95% доверительному интервалу.
 
Пятнадцать лет использования такого подхода к обработке сцинтилляционных спектров в десятках организаций Москвы и Московской области доказала его высокую эффективность и надёжность.
 

Поглощение излучения в веществе пробы

 
При прохождении через вещество гамма кванты взаимодействуют с электронами и ядрами атомов вещества. При взаимодействии они либо поглощаются в результате фотоэлектрического эффекта, либо теряют часть своей энергии и изменяют направление в результате комптоновского рассеяния и эффекта рождения электрон – позитронных пар (есть ещё ряд эффектов, но наиболее существенны три указанных). При малых энергиях  гамма квантов (до сотен кэВ), наиболее значимым является фотоэффект. При энергиях квантов Eγ ≥1022 кэВ, возникает эффект рождения электрон-позитронных пар. Влияние комптоновского рассеяния преобладает в промежуточной энергетической области. Совместное действие эффектов характеризуется полным линейным коэффициентом ослабления гамма квантов в веществе
 
     
Где n -концентрация атомов вещества в единице объема, σ -полное эффективное сечение взаимодействия.
      
Здесь ρ - плотность среды , Na - число Авогадро , A -молярная масса атомов среды.
      
       - сечения фотоэффекта, комптоновского рассеяния и эффекта образования пар.
В энергетическом диапазоне от 30 кэВ до 3000 кэВ
      ,       ,   
Z – заряд ядер атомов вещества счётного образца,  Eγ- энергия гамма излучения.
 
 
Рис. 4. Удельное эффективное сечение реакций взаимодействия гамма квантов с глиной (Z = 7.6).
 
 На рис. 4 представлены удельные эффективные сечения трёх основных реакций взаимодействия гамма квантов с веществом (кв.см/г) для счётного образца почвы (глины). Для энергий более 50 кэВ сечение фотопоглощения на порядок меньше сечения комптоновского рассеяния, сечение эффекта образования пар в диапазоне до 4000 кэВ также более чем на порядок слабее комптоновского рассеяния. Практически для основной массы исследуемых в Санэпиднадзоре веществ (продукты питания, почва и т.п.) единственным значимым эффектом является комптоновское рассеяние. Но встречаются пробы и с большими значениями Z эффективное:
 
Рис. 5. Удельное эффективное сечение реакций взаимодействия гамма квантов с вольфрамом (Z = 74).
 
На рис. 5 показаны аналогичные сечения реакций для образца вольфрама (в Московском центре гигиены и эпидемиологии пришлось исследовать образцы электродов из торированного вольфрама). В этом случае влияние фотоэффекта превышает влияние комптоновского рассеяния до энергии 500 кэВ.
 
Взаимодействие гамма излучения с веществом на пути от возникновения гамма кванта до регистрации в детекторе отражается в спектрах компонент рассчитываемых из эталонов при калибровке спектрометра. Но в комплексе счётный образец – детектор есть непостоянная часть – сам счётный образец. Поглощение излучения при прохождении через стенки измерительного сосуда, упаковку детектора и сам детектор при неизменной геометрии измерений остаётся постоянным, меняется поглощение излучения в веществе счётного образца. Это поглощение зависит от плотности, эффективной атомной массы и эффективного заряда атомов (Z эффективное) счётного образца.
 
Число гамма квантов на выходе из счётного образца, которые пролетят без взаимодействия с веществом пробы, выражается формулой:
                (1)
       N0– число гамма квантов для вакуума, x – эффективная толщина счётного образца.
 
Формула (1) вполне годится для расчёта поглощения гамма квантов в веществе счётного образца для метода обработки спектров связанных с выделением пиков (ППД спектрометрия), поскольку при любом взаимодействии гамма кванта с веществом этот гамма квант уже не зарегистрируется в пике полного поглощения – он будет потерян для обработки связанной с разложением спектра пиков (небольшое исключение составляет эффект образования пар, но этот эффект играет роль не в счётном образце, а в самом детекторе).
 
В сцинтилляционной спектрометрии (при обработке без выделения пиков) формула (1) в общем случае не работает – гамма кванты, испытавшие комптоновское рассеяние, могут не зарегистрироваться в детекторе, а могут и зарегистрироваться в спектре комптоновского рассеяния. Рассчитать, где в спектре зарегистрируется гамма квант (если зарегистрируется) после комптоновского рассеяния, очень непросто. В последнее время появилось много работ использующих для такого расчёта процедуру, основанную на разработанном в Los Alamos National Laboratory методе MСNP (Monte Carlo N-Particles transport code). К сожалению трудности связанные с корректным учётом комптоновского рассеяния остаются и при использовании этого метода – для учёта всех эффектов связанных с рассеянием гамма квантов требуется детальное описание не только системы счётный образец – детектор, но и защиты вокруг образца с детектором. Пока что мне не встречались работы убедительно доказывающие, что точность расчётов эффектов связанных с рассеянием гамма квантов с использование данного метода меньше, скажем, 10 %. Поэтому естественным выходом является перерасчёт (интерполяция или экстраполяция) спектров компонент на заданную плотность (комптоновское рассеяние пропорционально плотности электронов, а плотность электронов пропорциональна плотности счётного образца) по спектрам компонент, рассчитанным из эталонов различной плотности, приведённых к заданному эффективному заряду атомов и эффективной атомной массе. А перерасчёт на заданный эффективный заряд и эффективную атомную массу, то есть учёт фотопоглощения производить по формуле поглощения излучения в веществе (1). Именно такая процедура и реализована в программном пакете «Спектр».
 
Оценка плотности при измерениях в фиксированной геометрии проблем не вызывает – достаточно взвесить счётный образец, сложнее с Z эффективное и эффективной атомной массой пробы. Нередко оператор не знает химического состава пробы.
 
В пакете «СПЕКТР» предусмотрено несколько вариантов решения проблемы:
  • 1. Установка физических констант оператором.

  • 2. Определение Z эффективное по спектру пробы.

  • 3. Определение Z эффективное по спектру пробы и априорной информации об активности одного или нескольких нуклидов пробы.

 
Установка физических констант оператором
База данных по физическим константам (химические формулы вещества, эффективный атомный номер и эффективная атомная масса) в пакете «СПЕКТР» объединена с базой данных по нормативам (СанПиН, НРБ и т.п.). Оператор при запуске измерений должен выбрать из списка или создать тип пробы с требуемыми значениями физических констант. Форма процедуры создания – редактирования данных по физическим константам и нормативам выглядит следующим образом:
 
 
Norms 
Рис 6. Форма создания и редактирования констант и нормативов.
 
Оператор может ввести формулу вещества пробы для трёх вариантов счётных образцов. Задание нулевого значения эффективного атомного номера значит, что эти данные оператору неизвестны. В этом случае программа может попытаться определить требуемые данные непосредственно по измеренному спектру
 
Определение Z эффективное по спектру пробы
 В том случае если спектр пробы содержит нуклиды, имеющие гамма линии и в мягкой и в жёсткой области спектра, можно рассчитать Z эффективное непосредственно по спектру пробы, поскольку зависимость поглощения гамма излучения от энергии известна. Точность оценки будет тем выше, чем больше разница в активности нуклида, рассчитанная по мягким и жёстким линиям. Для включения процедуры расчёта Z эффективное по спектру пробы в процесс обработки спектра достаточно установить галочку «Варьирование Z эффективное пробы» в параметрах шаблона обработки.
 
 
 
Рис 7. Форма создания и редактирования параметров шаблонов декомпозиции.
 
На следующих двух рисунках представлены результаты обработки спектра электродов из торированного вольфрама измеренные на сцинтилляционном детекторе в Московском центре гигиены и эпидемиологии. Торий искусственного происхождения находится, как правило, в неравновесном состоянии, поэтому спектр раскладывается на две цепочки тория: Th-232 – Ac-228 и Th-228 – Tl-208.  
 
 
Рис 8. Результаты обработки спектра с варьированием Z эффективное.
 
Сравните результаты с обработкой того же спектра, но без варьирования Z эффективное (на следующем рисунке).
 
 
Рис 9. Результаты обработки спектра без учёта Z эффективное.
 
Значения активности Th-228 при обработке спектра без учёта Z эффективное существенно занижено, вследствие того, что не было учтен эффект фотопоглощения при большом Z эффективное пробы.
 
Определение Z эффективное по спектру пробы и априорной информации
В случае если проба не содержит нуклидов, имеющих линии и в мягкой (менее 100 кэВ), где значим процесс фотопоглощения, и в жёсткой (более 500 кэВ) области энергий, где он влияет слабо, автоматическое определение Z эффективное не работает, поскольку итоговая система уравнений становится вырожденной. В этом случае процедура решения системы уравнений исключает Z эффективное из списка неизвестных. Для того чтобы и в данном случае итоговая система уравнений не была вырожденной, можно дополнить её априорной информацией об активностях нуклидов имеющих линии в мягкой области энергий. Для этого потребуется либо химическое выделение требуемого нуклида и оценка его активности другими методами, либо добавление в счётный образец нуклида известной активности отсутствующего в пробе. В пакете «СПЕКТР» предусмотрена процедура учёта априорной информации об активностях нуклидов пробы либо вводимая оператором вручную, либо из файлов результатов обработки или архива измерений. Процедура решения итоговой системы уравнений использует помехоустойчивый вариант Байесовского оценивания.
 
 
8 октября 2011 г.
Дрёмин Геннадий Иванович.