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

Определение выходов ядерных реакций на основе анализа

цепочек распадов

С. С. Белышев1,0, К. А. Стопани2,0, С. Ю. Трощиев1, A.C. Курилик1, A.A. Кузнецов2

Московский государственный университет имени М. В. Ломоносова, 1 физический факультет, кафедра общей ядерной физики; 2 Научно-исследовательский институт ядерной физики имени Д. В. Скобельцына. Россия, И9991, Москва, Ленинские горы, д. /, стр. 2.

E-mail: " hatta@depni.sinp.msu.ru, ь belyshev@depni.sinp.msu.ru

Статья поступила 18.02.2011, подписана в печать 22.03.2011

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

Ключевые слова: ядерные реакции, активационный анализ, гамма-спектрометрия, линейные статистические модели.

УДК: 539.1.074.55. PACS: 25.20.^х.

Введение

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

к.

У (Ее) = ^П

cT(E)v(E,Ee)dE,

(1)

где Ее — энергия электронов, р(Е,Ее) — спектр тормозных фотонов, 1е — ток пучка, а п — число ядер в облучаемой части мишени. Для нахождения сечения реакции а(Е) необходимо провести несколько экспериментов с различной энергией электронов, после чего а(Е) определяется путем решения обратной задачи.

Во многих случаях непосредственное измерение выходов реакций путем регистрации ее продуктов невозможно (например, если требуется разделить каналы с вылетом различного числа нейтронов). Удобным методом измерений, позволяющим преодолеть данное затруднение, является метод наведенной активности [1]. В частности, данный метод является одним из основных при изучении многочастичных фотоядерных реакций (см., например, [2] и др. подобные работы). При использовании методики наведенной активности эксперимент состоит из двух этапов: на первом этапе исследуемый образец облучается в пучке тормозного излучения, а на втором переносится на детектор 7-квантов для измерения серии спектров остаточной активности. По пикам в измеренных спектрах идентифицируются распады нестабильных продуктов реакций и определяются их выходы.

1. Линейная модель площади пика в серии спектров

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

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

drij ~~dt

= yife(l) + Ц\]П] - XjO^ ,

Распад

i=l,...,N, (2)

Накопление

Рис. /. Пример цепочки распадов

'См. решение в частном случае, например, в [3].

где щ — количество ядер данного типа (или в данном изомерном состоянии), г/; — выход реакции образования г-го ядра, А/ — постоянные распада, Ц — вероятность распада /-го ядра в 1-е, /е(0 — ток ускорителя (если Тгш\ (время облучения), /е(£) = 0). Данная система может быть решена методом вариации постоянной, если переупорядочить уравнения таким образом, что в правой части все /<г, т.е. Ц= О при / > г. Это возможно всегда, так как в противном случае цепочка распадов была бы циклической. Тогда в первом уравнении системы все к/у = 0, и оно имеет решение

п\(0 = е

-А^

У\1е{т)е ]Тйт

(3)

,0

= е

У{1е(т)е 'т йт + щ

/=1 / (4)

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

N

/=1

(5)

где коэффициенты проще всего определяются,

если приравнять нулю все я? и все у, за исключением У! = 1. Коэффициенты сг(0 можно определить, положив все у-] равными нулю. Этим приемом можно воспользоваться как в численном расчете, так и в аналитическом.

Если ввести интегральный оператор

4/(0 = е'

еА"7(т) йт,

то сг-(0 можно представить в виде суммы

¿-1 с/

<*(') = ££ Я/*'

(6)

/=1 к=1

где для каждой последовательности / чисел

г .. > 7 > £ > £ > 1

= кГХМ%1 • • • ^Аг/г (я^^) . (7)

Аналогичная формула справедлива для апри этом выражение в скобках заменяется на уе1е1е(0-

г— 1

Общее число слагаемых в сг(0 равно 1 + С/ = 2г.

/=1

Отметим, что сг(0 равны нулю, если начальное количество ядер данного типа и тех, распады которых приводят к его образованию, равно нулю. Это так, если в состав облучаемой мишени не входят нестабильные ядра.

Измеряя серии из М спектров, определяют число распадов нестабильных ядер за промежутки времени Ы. Эта величина случайна и подчиняется закону распределения Пуассона РЙ(Х)> параметр которого ¡л, в спектре номер а равен

N

Ма ~ £ ЩаУ] + С" /=1

(8)

где

t¡,+At

tn+At

А1 —

Х{а1:(т)ёт, С1а =

Сг(т) йт.

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

Сделаем важное замечание. Несмотря на то что спектр 7-квантов, испускаемых при распадах различных ядер, обычно состоит из множества пиков и все они могут быть использованы для определения выходов, на практике использование данных о площадях всех пиков распада в рамках одной статистической модели затруднительно, хотя их энергии и относительные интенсивности известны. Дело в том, что эффективность детектора сложным образом зависит от энергии фотонов и от геометрии взаимного расположения детектора и образца. Существенное влияние на нее также оказывают эффекты истинного сложения пиков и самопоглощения фотонов в образце [4]. Вследствие этого соотношение между активностью образца и площадью пика в спектре известно лишь с точностью до коэффициента, не говоря уже о соотношении между площадями различных пиков одного и того же распада. Поэтому при обработке проще и надежнее определять набор выходов из анализа площадей разных пиков по отдельности. Результаты обработки обычно приводятся в виде, отнормированном на максимальный выход, поэтому при обработке пиков по отдельности неизвестный коэффициент сокращается. После анализа всех пиков будет получено несколько наборов выходов, которые можно рассматривать как несколько независимых измерений.

С учетом вышесказанного мы будем предполагать, что производится анализ некоторого пика распада г-го ядра и не будем в дальнейшем употреблять индекс г для обозначения г-го элемента цепочки. Из уравнения (8) следует, что числу распадов данного состояния за периоды измерения серии спектров соответствует линейная модель измерений, в которой математическое ожидание вектора числа распадов х равно

Ех Ау-с. (9)

где у — вектор неизвестных параметров г/;, а е — вектор случайных ошибок с нулевым математическим ожиданием и дисперсией, равной Ау. (В действительности, если то из вектора х необходимо вычесть вектор постоянных С, так что его распределение на самом деле Рца(ха — С). Это обстоятельство, однако, не меняет сути решения.)

Традиционным методом решения данной задачи в случае разных значений дисперсии вектора ошибок е является взвешенный метод наименьших квадратов [6], который позволяет найти несмещенную линейную оценку минимальной дисперсии у*:

у* = (А'П-1А)-1А'П

22 ВМУ. Физика. Астрономия. № 4

где О = diag(/ла) = diag(i?e) — диагональная вследствие независимости измерений ковариационная матрица ошибок. Покажем, что в случае распределения Пуассона данная оценка является наилучшей не только в классе линейных оценок. Найдем логарифм функции правдоподобия:

м

ln£ = In ( Ц

е '-'«//.f«

м

хп\

= + In /¿а ~ Щха\)).

\а=1 / а=1

Его производная по параметру г/; равна

дУ1 ¿V " Ма V

Найдем значения элементов матрицы информации Фишера:

(11)

(12)

hk(y)=E . „ „ \ду1 дук

dl_ dl

Vi

/мм , ч\

\а=1 /3=1

м м

vvfji Ak-Ai I E(xaxß) л} лй Л

a=lß=l

(13)

Для двух независимых Пуассоновских переменных ха и Xß

Е(х Xß)=[ C0V( A"" ' + E)CaExß = VaVß' если а i I3' \ох..-(Ех..)- /'..-/С- если a = ß.

(14)

Таким образом,

М А' Ак

!k ^ ßa

(15)

а=1

что позволяет избежать случаев, когда ха = 0. Кроме того, когда имеет место распределение Пуассона, оценка (1 /¡л)* = \/(х + \) имеет меньшее смещение, чем оценка \/тах(х,\) (смещение предлагаемой оценки меньше 10% при х>5, что легко показать численно). Метод, в котором данная оценка дисперсии используется для получения оценки параметров

y-FGLS={ATn-i:A)-lATn-l1x, (17)

называется доступным методом наименьших квадратов (feasible generalized least squares (описание метода и ссылки на литературу см., например, в [8]). Добавим, что именно в этой форме он часто применяется для анализа серий спектров в гамма-активационных экспериментах.

Другим способом получения оценок в обобщенных линейных моделях является метод наименьших квадратов с итерационным перевзвешиванием (iteratively reweighted least squares), в котором полученная на каждом шаге оценка у* используется для уточнения оценки О^1* = й\щ{Ау*п)^1, т.е. следующая оценка У:. 11 равна

у*п+1 = (Аг А[щ{Ау*пу{АГ1АГ 6.mg(Ay*ny{X. (18)

В [7] описано использование данного метода для получения оценок максимального правдоподобия в моделях, принадлежащих к экспоненциальному семейству распределений. Покажем, как данный метод может быть использован в нашем случае. В качестве начальной оценки можно воспользоваться, например, несмещенной оценкой обычного МНК, в который ковариационная матрица не входит. На каждой итерации происходит уточнение весов, и элементы вектора х с большей величиной получают больший вес. Схождение данного процесса соответствует следующему предельному условию:

■1 /П^1 лТ

у*00 = (А< diagiAy^y'A) 1ЛЧ1аё(Л14) х, (19)

или

ATdiag(Ayl0ylAyl0 = AT-l=ATdiag(Ayl0ylx, (20)

или в матричном виде I = АтО,^1А. В то же время дисперсия оценки по методу МНК Оу* = ,

что совпадает с нижней границей неравенства Крамера-Рао Оу* и данная оценка, следовательно, эффективна.

2. Оценки выходов в случае неизвестной ковариационной матрицы

Описанный в предыдущем разделе метод непригоден для обработки результатов эксперимента, так как значения матрицы О неизвестны. Это затруднение характерно и для других задач линейной регрессии с неодинаковыми дисперсиями. Обычно оно преодолевается применением метода наибольшего правдоподобия или использованием некоторой оценки матрицы О^1*, подставляемой в (10).

В качестве оценки О^1* можно использовать матрицу

= diag(l/(л;a + 1)) (вместо А\&%(\/ха)), (16)

или

м м А

(21)

Таким образом, в случае существования предела lim у* у* его значение совпадает с оценкой мак-

п—»оо

симального правдоподобия (12). Кроме того, в пределе больших ха значения у^ и оценки максимального правдоподобия совпадают с оценкой взвешенного МНК, что соответствует хорошо известному результату: в случае нормального распределения ошибок в линейной модели оценка максимального правдоподобия равна оценке взвешенного МНК и эффективна [5].

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

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

Бу* = (ЛгГГ>Л)^ = (Лг(0^1)1/2(0^1)1/2Л) = (ВТВ),

(22)

где (О ')1/2(0 ')1/2 — разложение Холецкого матрицы О^1, В — нормализованная матрица, соответствующая линейной модели с единичными дисперсиями измерений, соответствующий ей вектор измерений — (О^1)172*- Сложению части спектров серии соответствует новая матрица С, строки которой соответствуют суммам строк матрицы В с последующей нормировкой. Рассмотрим ситуацию, когда суммируются последние гп строк матрицы В. Это допущение является позволительным, так как суммируемые строки перестановкой могут быть перемещены в нижнюю часть матрицы В, что не скажется на результате. Введем следующие обозначения:

В =

Ву

Вт

с =

(23)

где Вт — последние т строк матрицы В, В\ — верхняя часть этой матрицы, Ст — строка, состоящая из сумм столбцов матрицы Вт, нормированной на 11\[ш. Дисперсия оценок выходов в этом случае будет определяться обратной матрицей (С1"С)-1, не зависящей от перестановки строк в матрице В. Рассмотрим элементы матриц В1В и СТС:

М—т— 1

М

{ВТВ)а= £ вк(Вк1+ £ Вк1Вщ,

к=М-т

к= 1 М—т— 1

{СТС)а= £

к= 1

1

т

м

м

Е М Е Д

Кк=М-т ) \к=М—т

(24)

(25)

(26) (27)

Можно записать в более краткой форме: ВТВ = В\В{+ВТтВт, СТС = В[В\ + с1ст.

Рассмотрим разность матриц В В — С С = ВтВт — СтСт = = ВТт1тВт - ВТт11ТВт = ВТ - вт, (28)

где 1т — единичная матрица тх т, 1 — единичный вектор. Известно, что для любого вектора а матрица 1т - ааТ > 0, если ата [6]. В нашем случае ата = г 1 1. т.е. данное условие выполняется1. Далее, если некоторая матрица Р ^ 0, то С}ТРС} > 0 для всех матриц С). Из этого следует, что ВтВ - СТС ^ О, или ВтВ > СТС. В соответствии с [9] это значит, что план эксперимента, задаваемый матрицей В, более оптимален, чем план, задаваемый матрицей С. Единственным случаем, когда увеличение дисперсии оценок не происходит, является суммирование измерений, соответствующих одинаковым строкам В. Это понятно и на качественном уровне: одинаково распре-

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

Приведем результаты моделирования, демонстрирующие применение описанных методов. Рассматривалась цепочка, состоящая из четырех элементов с постоянными распада 2.0, 0.1, 0.7 и 0.0. Каждый элемент с одинаковой вероятностью распадается во все последующие. На основе этих данных была рассчитана матрица А, соответствующая 100 измерениям спектров с начального момента до 1= 100. Четвертый столбец матрицы равен нулю, поэтому был отброшен. Результаты определения выходов разными методами приведены в таблице.

Результаты моделирования позволяют, в частности, сделать вывод о том, что обработка данных методом ЕОЬБ (доступный МНК) дает несколько худшие результаты, чем другие методы.

3. Особенности анализа реальных данных

Для того чтобы применять описанные оценки на практике, требуется сделать несколько уточнений. Первое связано с тем, что в реальности измеряется не число распадов, а площадь пика, имеющего форму функции Гаусса. Обычно пик содержится в нескольких соседних каналах спектра, каждый из которых представляет собой независимое измерение. Поэтому число уравнений в линейной модели (9) увеличивается в соответствующее число раз, т. е. уравнение номер а заменяется на Ь уравнений:

АафУ = хаф, I 1.....Ь. А1а[}=А1а{с1{Е[м)^с1{Е())),

(29)

где ¿(£г) — ¿(Е\) — интеграл функции Гаусса от Е\ до £2 • Эта замена тривиальна, поскольку, не говоря уже об энергии пика, ширина пика при данной энергии хорошо известна пользователю детектора и с большой точностью стабильна при разных скоростях счета.

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

Что касается искусственной части фона, то ее можно учесть путем введения дополнительных параметров в линейную модель. На рис. 2 показан отклик детектора на поток 7-квантов с энергией 662 кэВ (стандартный образец 137 Сэ). Он состоит собственно из пика полного поглощения и области континуума. Площадь обеих областей пропорциональна интенсивности распадов, а форма в основном определяется особенностями детектора и хорошо описывается моделированием (см., например, [10, 11]). В реально измеренных спектрах

Запись А ^ 0 означает, что матрица А положительно полуопределенна, а А^ В значит, что А — В ^ 0.

Результаты определения выходов по сгенерированным данным*

i OLS а, % WLS a, % FGLS а, % IGLS а, % у/Щ*> %

I. Ух = 103, y2 = 104, ¿/3 = 2- 104. Площадь пика от 150 до 7000

1 1 003.39 74.12 1 004.99 64.88 1 036.43 64.94 1 004.98 64.88 64.43

2 9 998.49 5.10 9 996.91 4.05 9 930.63 4.05 9 996.91 4.05 4.04

3 19999.15 0.82 19 998.95 0.76 19991.91 0.76 19 998.95 0.76 0.76

II. yx = 105, y2 = 106, £/з = 2- 106. Площадь пика от 15000 до 700000

1 99 945.41 7.38 99 969.94 6.40 100001.57 6.40 99 969.95 6.40 6.44

2 1000045.76 0.51 1000019.78 0.40 999 953.51 0.40 1000019.77 0.40 0.40

3 2000013.12 0.08 2 000011.03 0.08 2 000003.96 0.08 2 000011.03 0.08 0.08

* Приведены средние значения и среднеквадратичные отклонения оценок, полученных разными методами. Число повторений: 10000. Последний столбец — оценка стандартного отклонения взвешенного МНК.

Число отсчетов 1е+06 г

100 000

10 000 :

1 000

300 400 Энергия, кэВ

700

Рис. 2. Отклик детектора на поток монохроматических 7-квантов

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

+ ЕафУЬ = Ха,/3,

(30)

где Е]ар — вклад /-го высокоинтенсивного распада в фон в канале спектра /3 в измерении номер а, а уь — выходы реакций, в которых рождаются эти ядра. Число элементов в у^, как уже было сказано, невелико, а коэффициенты можно определить,

например, с помощью метода Монте-Карло. Таким образом, системой (30) полностью описывается несколько каналов спектра в окрестности исследуемого пика, если наборы выходов у и уь не пересекаются. При условии достаточно полного знания формы спектра в области континуума метод решения, описанный в предыдущем

разделе, позволяет определить неизвестные выходы у при наличии фона.

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

Список литературы

1. Кузнецов Р.А. Активационный анализ. Изд. 2-е. М., 1974.

2. Алиев Р.А., Ермаков А.Н., Ишханов Б.С. и др. // Вестн. Моск. ун-та. Физ. Астрон. 2006. № 6. С. 55.

3. Bateman Н. // Proc. Cambridge Philos. Soc. 1910. 15. P. 423.

4. Gilmore G.R. Practical gamma-ray spectrometry. 2nd ed. John Wiley & Sons, 2008.

5. Боровков А.А. Математическая статистика. Новосибирск, 1997.

6. Rao C.R., Toutenburg H., Shalabh, Heumann C. Linear models and generalizations. Springer, 2008.

7. McCullagh P., Nelder J.A. Generalized linear models. Chapman and Hall, 1989.

8. Beck N., Katz J.N. // Amer. Polit. Sci. Rev. 89, N 3. P. 634.

9. Математическая теория планирования эксперимента /

Под ред. С.М. Ермакова. М., 1983. 10. Трощиев С.Ю. Ц Труды X межвузовской научной школы молодых специалистов «Концентрированные потоки энергии в космической технике, электронике, экологии и ме-

дицине» / Под ред. Б. С. Ишханова, Л. С. Новикова. М., 2009. С. 174.

11. Berlizov A.N., Magill J. An interactive Web accessible garnrna-spectrurn Simulator. http://www.nucleonica.net:81/ wiki/irnages/0/03/Nucleonica3.pdf.

Measuring nuclear reactions yields in a decay chain analysis based procedure

S.S. Belyshev1 ,b, K. A. Stopani2 fl, S.Yu. Troschiev1, A.S. Kurilik1, A. A. Kuznetsov2

1D. V. Skobeltsyn Institute of Nuclear Physics, M. V. Lomonosov Moscow State University, Moscow 119991, Russia.

2 Department of General Nuclear Physics, Faculty of Physics, M. V. Lomonosov Moscow State University, Moscow 119991, Russia.

E-mail: a hatta@depni.sinp.msu.ru, b belyshev@depni.sinp.msu.ru.

The problem of finding nuclear reactions yields in activation experiments is being analysed. It was shown that in the presence of complex decay chains the activity of gamma-spectrum lines of the nuclei can be described with a linear statistical model. Since the covariance matrix is unavailable the iteratively reweighted general least squares method is an appropriate means of obtaining a maximum-likelihood estimation. The relation between the variances of estimated yields and the duration of spectrum measurements is being discussed.

Keywords: nuclear reactions, activation analysis, gamma-ray spectrometry, linear statistical models.

PACS: 25.20.-x.

Received 18 February 2011.

English version: Moscow University Physics Bulletin 4(2011).

Сведения об авторах

1. Белышев Сергей Сергеевич — физик; тел.: (495) 939-25-58, e-mail: belyshev@depni.sinp.msu.ru.

2. Стопани Константин Александрович — мл. науч. сотр., тел.: (495) 939-25-58, e-mail: hatta@depni.sinp.msu.ru.

3. Трощиев Сергей Юрьевич — аспирант, тел.: (495) 939-56-35, e-mail: sergey.troschiev@googlemail.com.

4. Курилик Александр Сергеевич — ассистент; тел.: (495) 939-56-35, e-mail: kurilik@depni.sinp.msu.ru.

5. Кузнецов Александр Александрович — мл. науч. сотр.; тел.: (495) 939-25-58, e-mail: herero@mail.ru.