Спросить
Войти
УДК 519.245+51-7
ПРИМЕНЕНИЕ МЕТОДА МОНТЕ-КАРЛО ДЛЯ РАСЧЕТА НЕОПРЕДЕЛЕННОСТИ ДОЗЫ ОБЛУЧЕНИЯ ЩИТОВИДНОЙ ЖЕЛЕЗЫ
С.С. ТРЕТЬЯКЕВИЧ
Государственное учреждение образования «Белорусская медицинская академия последипломного образования»
(БелМАПО), г. Минск
Введение
С 1990 года на территории Беларуси наблюдается рост числа заболеваний щитовидной железы (ЩЖ) [1], обусловленный радиационным облучением ЩЖ вследствие аварии на Чернобыльской АЭС (ЧАЭС) [2]. Неопределенность величины радиационного риска заболеваний ЩЖ потребовала дополнительных исследований зависимости «доза - эффект». Одно из исследований проводится на территории Беларуси с 1996 года в рамках Белорусско-Американского научного Протокола по изучению рака и других заболеваний ЩЖ в Беларуси после аварии на Чернобыльской АЭС [3].
Основным источником формирования дозы облучения ЩЖ жителей Беларуси
после аварии на ЧАЭС являлось внутреннее облучение ЩЖ I [4]. Наиболее
достоверно доза облучения ЩЖ I может быть реконструирована на основе
прямых инструментальных измерений мощности экспозиционной дозы (МЭД) над
ЩЖ в период действия радиоактивного йода («йодный» период). По результатам
прямых измерений МЭД над ЩЖ жителей Беларуси, выполненных в мае-июне 1986
года, и типовым условиям жизнедеятельности были определены индивидуальные
дозы облучения ЩЖ I на 130 тысяч человек [5]. Используя предварительные данные оценки доз, в когорту Белорусско-Американского научного Протокола были отобраны примерно 12 тысяч человек, имевших на момент измерения возраст до 18 лет.
У большинства субъектов когорты МЭД над ЩЖ в течение «йодного» периода была измерена однократно. При оценке дозы облучения ЩЖ это обстоятельство
требует привлечения дополнительных сведений о поступлении I в организм субъекта ингаляционным и пероральным путем в «йодный» период (места пребывания субъекта, источники и количество потребления молока, молочных продуктов и зеленых листовых овощей, употребление препаратов стабильного йода).
Измерения МЭД над ЩЖ жителей Беларуси были выполнены неспециализированными приборами (ДП-5, СРП-68-01) и во многих случаях
непрофессиональным персоналом. В связи с этим оценка доли МЭД над ЩЖ от I, содержащегося в ЩЖ, характеризуется значительной неопределенностью. Данные о
поступлении I в организм субъекта ингаляционным и пероральным путем в «йодный» период, получаемые путем индивидуального опроса субъекта, также вносят свой вклад в погрешность оценки дозы из-за ненадежности ответов субъекта на вопросы индивидуального опроса, вызванной большим периодом времени, прошедшим после аварии на ЧАЭС. На результирующую неопределенность дозы
облучения ЩЖ I будет влиять и вариабельность параметров модели оценки дозы, обусловленная природной вариабельностью параметров или незнанием их подлинных значений.
Постановка задачи
При вычислении дозы облучения ЩЖ субъектов когорты, обусловленной I, используются два подхода. В работе [6] подробно описана модель реконструкции
дозы облучения ЩЖ I при отсутствии прямых измерений МЭД над ЩЖ путем радиоэкологического моделирования. Рассчитанная таким способом доза облучения ЩЖ 1311 далее обозначается термином «экологическая» доза. В работе [7]
представлена модель расчета дозы облучения ЩЖ I при наличии прямых измерений МЭД над ЩЖ. Доза, рассчитанная по данным прямых измерений МЭД над ЩЖ, далее обозначается термином «инструментальная» доза.
Большинство параметров модели реконструкции дозы облучения ЩЖ I характеризуются вариабельностью, выражаемой распределениями разного типа [6],
В сумме модель реконструкции дозы внутреннего облучения ЩЖ I содержит 42 параметра, значения которых являются вариабельными и характеризуются 4-мя типами распределений: 13 параметров описываются равномерным
распределением, 10 - треугольным (Симпсона), 1 - нормальным, 18 логнормальным [6], [7]. Наличие распределений разного типа не дает возможности
получить явное выражение для расчета неопределенности дозы облучения ЩЖ I. Единственным способом вычисления неопределенности является использование метода статистического моделирования (метод Монте-Карло).
В используемых моделях реконструкции доз облучения содержание 131! в ЩЖ определяется через систему из четырех камер, описываемых дифференциальными уравнениями первого порядка:
= ¿(г) ‘ ^1(ґ)
= Л(02 • У 2 (0 й (1)
= У3( у^ у 2, г)-*< 3 • У3(г)
йУ7^) = Л(У^У2,Уз,г)4 • У4(гX йг
где уі - содержание 131! в камере і в момент времени г, кБк; /і - поступление 131! в камеру і в момент времени г, кБк-сут"1; Хі -постоянная выведения 131! из камеры і, сут-1.
В применяемых моделях камерами являются: 1 - пастбищная трава; 2 -пастбищная почва; 3 - молоко частной коровы; 4 - ЩЖ.
Доза облучения ЩЖ I определяется как интеграл содержания I в ЩЖ и выражается следующей формулой:
D = DF • | у4(г)йг,
где D - поглощенная доза облучения ЩЖ I за период времени АГ, мГр; DF - доза
на единицу распада I в ЩЖ, мГр.
Таким образом, система дифференциальных уравнений и интегралы делают невозможным применение для вычисления неопределенности такого математического пакета, как Crystal Ball, обрабатывающего только алгебраические выражения. Применение математических пакетов MathCad или MatLab требует больших затрат времени по переводу большого количества данных индивидуального опроса в формат, используемый MathCad или MatLab.
Наилучшим выходом в такой ситуации является разработка прикладного
программного обеспечения для вычисления индивидуальных доз облучения ЩЖ I и их неопределенностей. Основной задачей при разработке такого программного обеспечения является корректная реализация метода Монте-Карло. В настоящей работе для расчета доз и неопределенностей использовались возможности, предоставляемые системой управления базами данных (СУБД) Microsoft Access 8.0 и встроенным в эту СУБД языком программирования Visual Basic for Application (VBA) 8.0. Выбор языка программирования был обусловлен хранением всех исходных данных в формате СУБД Microsoft Access 8.0.
Реализация метода Монте-Карло
Сущность метода Монте-Карло заключается в многократном вычислении выходных параметров системы под воздействием входных параметров с заданными законами распределения. Для каждого из параметров в соответствии с функцией плотности вероятности параметра генерируется случайное значение. По множеству полученных случайным образом значений параметров рассчитываются выходные параметры (в данном случае выходной параметр только один — доза облучения
ЩЖ I). Повторив подобную процедуру N раз, получим распределение N значений дозы облучения ЩЖ 131I, из которого рассчитывается неопределенность оценки дозы.
Основой реализации метода Монте-Карло является генерация последовательности случайных чисел с равномерным распределением в интервале {0, 1}. В настоящей работе для получения такой последовательности используются два способа генерации случайных чисел. Один способ основан на встроенном в VBA 8.0 генераторе случайных чисел, который инициализируется значениями, возвращаемыми системным таймером. Другой способ получения последовательности случайных чисел основан на использовании генератора псевдослучайных чисел, в качестве которого используется стандартный минимальный генератор с периодом (231-2) [8]. Он генерирует квазиравномерно распределенные целые числа в диапазоне {0, 231-3}. Реализация данного генератора в форме VBA-функции представлена в Приложении 1. Переход к последовательности случайных чисел, равномерно распределенных в интервале {0, 1}, производится путем деления числа, возвращаемого генератором случайных чисел, на значение верхней границы диапазона возвращаемых чисел (231-3).
Получение последовательностей случайных чисел, распределенных по произвольному закону, осуществляется путем их формирования из исходной последовательности случайных чисел, равномерно распределенных в интервале {0, 1}, с помощью различных формул преобразования.
Переход от последовательности равномерно распределенных в интервале {0, 1} случайных чисел Vi к последовательности случайных величин x, c равномерным законом распределения на отрезке {a, b} выполняется по формуле [9]:
х = а+(Ь - а) •у, , (3)
где а - нижняя граница отрезка; Ь - верхняя граница отрезка.
Преобразование последовательности равномерно распределенных в интервале {0, 1} случайных чисел Vі к последовательности случайных величин х, с
треугольным законом распределения на отрезке {а, Ь} и модой в точке т, не совпадающей с точкой (а+Ь)/2, производится по формуле [8]:
а + уІV, ■ (т - а) • (Ь - а), х <
_______________________________________________________ Ь - а
І т - а
Ь - (т - Ь) • ((т - Ь) + (Ь - а) • (V, - ----------------------------)), х > —И
Ь - а Ь - а
Переход от последовательности равномерно распределенных в интервале {0, 1} случайных чисел Vi к последовательности случайных величин xi с нормальным законом распределения осуществляется по методу Марсальи-Брея (Marsaglia-Bray). Реализация этого метода в виде "УЪА-функции, выполненная в соответствии с работой [8], представлена в Приложении 2.
Преобразование последовательности равномерно распределенных в интервале {0, 1} случайных чисел Vi к последовательности случайных величин xi с
логарифмически нормальным законом распределения производится на основе метода Марсальи-Брея. Реализация данного перехода в виде "УВА-функции, выполненная согласно [8], показана в Приложении 3.
Вычисление неопределенности дозы облучения ЩЖ I, реализованное в программном коде, производится в следующей последовательности:
каждого из параметров модели реконструкции дозы облучения ЩЖ I, значения которых являются вариабельными:
- генерация опорного числа для генератора псевдослучайных чисел (для каждого параметра используется другое опорное число);
- передача опорного числа в генератор псевдослучайных чисел;
- генерация последовательности равномерно распределенных в интервале {0, 1} псевдослучайных чисел мощностью N
- преобразование полученной последовательности равномерно распределенных в интервале {0, 1} псевдослучайных чисел мощностью N в последовательность случайных чисел с законом распределения параметра.
интерполирование и интерполирование кубическими сплайнами реализовано в программном коде с помощью численных методов. Вычисление интегралов производится по общей формуле Симпсона [10]. Решение системы дифференциальных уравнений 1-го порядка осуществляется по методу Рунге-Кутта четвертого порядка [9]. Интерполирование выполняется методом интерполяции кубическими сплайнами [9] и методом кусочно-линейной интерполяции [10]. В компьютерной программе предусмотрен розыгрыш 1000, 10000 и 100000 историй
для вычисления результирующей неопределенности дозы облучения ЩЖ I.
Результаты
В базе данных индивидуального опроса хранятся результаты первичного анкетного опроса для 11767 субъектов когорты. Для каждого из этих субъектов просчитано 1000 историй и получено 1000 значений дозы облучения ЩЖ 131Г Время расчета составило 2,5 суток на персональном компьютере следующей конфигурации: процессор - АШІоп ХР 1700+; оперативная память — 256 МВ; жесткий диск — 20 GB.
Выборки из 1000 значений дозы облучения ЩЖ I описываются распределением логарифмически нормального типа. Это позволяет выражать неопределенность оценки дозы облучения ЩЖ 131! через стандартное геометрическое отклонение (СГО). Распределение субъектов когорты в зависимости от СГО, определенного по выборке из 1000 значений дозы облучения ЩЖ 131^ представлено на рис. 1 и 2. На рис. 1 представлено распределение 11165 субъектов когорты с «инструментальной» дозой, на рис. 2 - 602 субъектов когорты с «экологической» дозой.
Рис. 1. Распределение СГО для субъектов когорты с «инструментальной» дозой
Значения СГО, полученные для субъектов когорты с рассчитанной «инструментальной» дозой облучения ЩЖ I, находятся в интервале 1,66 - 4,20. Медиана выборки, составленной из значений СГО для 11165 субъектов когорты с «инструментальной» дозой, равна 2,0.
Ё 70 ш
Рис. 2. Распределение СГО для субъектов когорты с «экологической» дозой Значения СГО, полученные для субъектов когорты с вычисленной
«экологической» дозой облучения ЩЖ I, находятся в интервале 1,85-3,46. Медиана выборки, составленной из значений СГО для 602 субъектов когорты с «экологической» дозой, равна 2,3.
Полученные результаты сходятся с данными работы [11], в которой с помощью математического пакета Crystal Ball были оценены неопределенности доз облучения
I ЩЖ жителей Беларуси, рассчитанных путем радиоэкологического моделирования. Анализ, выполненный по 20 параметрам модели, показал, что численные значения фактора неопределенности в единицах СГО составляют 2,3-2,4 для средних по населенному пункту и 2,6-3,1 для индивидуальных доз облучения ЩЖ1311[11].
Заключение
Средствами языка программирования VBA 8.0 разработан инструмент для расчета
неопределенности доз облучения ЩЖ I, позволяющий за приемлемое время рассчитать дозы облучения и их неопределенности для большого количества субъектов.
Трехмодовый вид распределения СГО, показанный на рис. 1, обусловлен тем, что неопределенность прямых измерений МЭД над ЩЖ субъектов когорты описывается только 3 значениями СГО: 1,4, 1,8, 2,3. Рисунок 1 показывает, что вариабельность данного параметра является определяющей при расчете
неопределенности «инструментальной» дозы облучения ЩЖ 131I и для ряда субъектов может быть завышена.
Assembly, with Scientific Annexes. Annex J: Exposure and effects of the Chernobyl accident. - New York: United Nations, 2000. - P. 451-566.
Приложение 1
VBA-функция реализации стандартного минимального генератора псевдослучайных чисел
Function GetRN (IngSeed As Long) As Long Const IngA As Long = 16807 Dim IngM As Long Dim IngQ As Long Dim IngR As Long IngM = 2 A 31 - 1 IngQ = IngM \\ IngA IngR = IngM Mod IngA
IngSeed = IngA * (IngSeedMod IngQ) - IngR * (IngSeed \\ IngQ)
If IngSeed < = 0 Then
IngSeed = IngSeed + IngM End If
GetRN = IngSeed - 1 End Function
Приложение 2
VBA-функция реализации перехода к последовательности нормально распределенных чисел
Function GetNormalRN (sngMAs Single, sngSD As Single) As Double Dim R1, R2, R3, S As Double Randomize Do
R1 = 2# * Rnd - 1#
R2 = 2# * Rnd- 1#
R3 = R1 * R1 + R2 * R2 Loop While R3 < 1#
S = Sqr (Abs(2# * Log(R3) /R3))
GetNormalRN = sngM + sngSD * S * R1 End Function
Приложение 3
VBA-функция реализации перехода к последовательности логарифмически нормально распределенных чисел
Function GetLogNormalRN (sngGM As Single, sngGSD As Single) As Double Dim R1, R2, R3, S As Double Dim A, B, C, D As Double Randomize Do
R1 = 2# * Rnd- 1#
R2 = 2# * Rnd- 1#
R3 = R1 * R1 + R2 * R2 Loop While R3 < 1#
A = (Log(sngGSD)) A 2 B = Exp(Log(sngGM) + A / 2)
C=-A/2 D = -A * 2
S = Sqr(Abs(D * Log(R3) /R3))
GetLognormalRN = B * Exp(C + S * R1)
End Function
Получено 16.06.2005 г.
