CC BY
59
ПРИКЛАДНАЯ МАТЕМАТИКА И МЕТОДЫ МАТЕМАТИЧЕСКОГО МОДЕЛИРОВАНИЯ
УДК 621.039.9; 621.396
А. А. Г р е ш и л о в, А. Л. Лебедев, П. А. Плохута
ГАЗООБРАЗНЫЕ ПРОДУКТЫ ДЕЛЕНИЯ И СЕЙСМИКА КАК ИДЕНТИФИКАТОРЫ ЯДЕРНЫХ ВЗРЫВОВ
Рассмотрена задача идентификации ядерных взрывов по газообразным продуктам деления, изотопам криптона и ксенона, и по сейсмическим измерениям. По измеренным в атмосфере активностям изотопов благородных газов оценивается наличие или отсутствие подземного ядерного взрыва (в идеальном случае — вклады разных источников радиоактивности в суммарную измеренную активность изотопа), а на основе данных сейсмических станций (задача пеленгации) — географическое положение возможной зоны проведения ядерных испытаний. Приводятся алгоритмы получения точечных и интервальных оценок решения.
Ключевые слова: идентификация, регуляризация, математическое программирование, изитоп, сейсмика, конфлюэнтный анализ.
1. Идентификация по измеренным активностям изотопов криптона и ксенона. Применяемые методы идентификации. В процессе ядерного взрыва происходит мгновенное деление ядер 235и, 238и, 239Ри. Под мгновенным делением понимается одномоментное деление тяжелых ядер (деление в процессе взрыва), в результате которого возникают осколки деления, распадающиеся с течением времени по цепочкам радиоактивных превращений (получаются продукты деления). Количество осколков деления определяется их независимыми выходами, количество продуктов деления определяется кумулятивными или полными выходами данной цепочки. При ядерном взрыве в атмосферу могут попадать изотопы Кг и Хе. Задача идентификации по измерениям активностей изотопов этих инертных газов должна решаться с учетом помех — возможных выбросов аналогичных изотопов в атмосферу при функционировании ядерных реакторов, импульсных источников нейтронов на основе делящихся тяжелых ядер, при переработке отработанного ядерного горючего. Во всех этих случаях в атмосферу попадают радиоактивные изотопы благородных (инертных) газов — продукты деления изотопов тяжелых ядер (урана и плутония). В итоге в пробах, отобранных в атмосфере, суммарная активность изотопов Кг и Хе может быть обусловлена наличием нескольких разных источников.
Изотопы Кг более чувствительны к различным видам деления (выход Кг зависит и от делящегося материала, и от энергии делящих нейтронов). Изотопы Хе мало чувствительны к различным видам деления и потому идентификация по ним представляет собой более трудную задачу.
Хотя с течением времени спустя 24 ч с момента мгновенного деления активность определяется в основном изотопами Хе (см. табл. 13 в работе [1]), но поскольку для мониторинга радиоактивной обстановки по всему земному шару создается сеть станций, то может случиться, что источник выброса окажется недалеко от регистратора и будут зарегистрированы изотопы и Кг, и Хе. Более детальные измерения позволяют определить концентрацию изотопов Кг и по прошествии 24 ч.
Будем рассматривать алгоритмы решения общей задачи, когда измерены активности изотопов Кг и Хе, а частная задача — идентификация ядерного взрыва по изотопам Хе — получается из общей за счет уменьшения числа переменных — зарегистрированных изотопов.
В СССР метод определения параметров ядерного взрыва по изотопам Кг и Хе был предложен А.А. Грешиловым и защищен авторскими свидетельствами (авт. свид. 43815 с приоритетом 4 марта 1967 г., авт. свид. 49143 с приоритетом 10 февраля 1968 г., авт. свид. 46128 с приоритетом 29 апреля 1968 г.), способ определения активности газов Кг и Хе в смеси газообразных продуктов деления изложен в изобретении [2]. В отобранных пробах воздушной среды определялась активность Кг и Хе, а затем рассчитывался относительный вклад различных видов деления в энерговыделение данного взрыва. Активность отобранной пробы определялась в течение недели, т.е. имитировалась ситуация с распространением радиоактивной струи. Эксперименты проводились в 1968 г. на Семипалатинском полигоне.
В работе [3] подробно описан математический алгоритм решения некорректной задачи оценки вкладов различных видов деления в общее энерговыделение.
Одним из предлагаемых в последнее время методов идентификации источников изотопов благородных радиоактивных газов является метод, в котором используются отношения активностей изотопов Хе: 131тХе, 133тХе, 133Хе и 135Хе [4].
Суть метода заключается в том, что для разных источников радиоактивных газов (ядерный взрыв, реакторный выброс, фоновая активность) независимо друг от друга строятся графики зависимости отношений активностей одних изотопов Хе от отношений активностей других изотопов Хе. Один из таких графиков приведен на рис. 1.
10 5 10° 105 ШтХе/131тХе
Рис. 1. Графики зависимости отношения активностей 135Xe/133Xe от 133mXe/131mXe для случая ядерного взрыва, реакторной и фоновой активностей
Поскольку между областями реакторной, фоновой активности и ядерных взрывов может существовать разделяющая линия (области не налагаются друг на друга), то предлагается [4] для измеренных активностей изотопов Хе строить эти отношения и по их принадлежности к той или иной зоне говорить о наличии (либо отсутствии) ядерного взрыва и реакторного выброса.
В реальных условиях в отобранной пробе могут быть изотопы Хе от различных источников, включая фоновую активность этих изотопов. Поэтому отношения активностей изотопов Хе будут определяться вкладом каждого источника: фон плюс ядерный взрыв могут идентифицировать ядерный реактор.
Кроме того, активность изотопов в отобранной пробе, а соответственно, и отношение активностей в значительной мере зависят от процессов сепарации изотопов Хе от предшественников по цепочкам радиоактивных превращений.
Следовательно, предложенный в [4] подход является идеализированным и в реальных условиях не применим по нескольким причинам:
1) если в атмосфере присутствуют радиоактивные изотопы Кг и Хе различного происхождения, то непонятно, какой вклад каждый источник радиоактивных изотопов внес в общую активность и, соответственно, в отношения активностей, а потому невозможно однозначно построить отношения для изотопов от разных источников;
2) поскольку графики строятся для отношения активностей, то теряется часть информации, так как относительные единицы являются "вторичным" продуктом по сравнению с абсолютными;
3) отношение активностей изотопов Кг и Хе в ядерном реакторе зависит от времени кампании, т.е. может происходить смещение линий из одной зоны в другую (например, реактор можно "продуть" раньше, отчего реакторная зона перемещается в зону ядерных взрывов);
4) в струе выходящих из места взрыва газов (в случае проведения подземного взрыва в штольне или скважине) можно провести смешивание изотопов, препятствующее определению настоящего источника активности (например, организатор ядерных испытаний может "добавить" в зону проведения ядерного взрыва определенное количество изотопов Хе того или иного вида и тем самым завуалировать источник выброса радиоактивных изотопов);
5) не учитывается время сепарации изотопов Кг и Хе от их предшественников по цепочкам радиоактивных превращений, что дает неправильные представления об изменения активности изотопов во времени и приводит к неправильным выводам.
На рис.2 приведен пример влияния времени кампании реактора на отношения активностей (на рис. 1 нанесены значения отношений активностей изотопов Хе в случае работы реактора с временем кампании от 0,01 года до 3 лет [5]). Зона реакторной активности (красные линии) переместилась в зону ядерных взрывов.
По указанным причинам для идентификации источников радиоактивности изотопов благородных газов необходимо исследовать только абсолютные значения измеренных активностей, учитывать время се-
10'
10'
10
-5
10
Время
Вред кампа реакт« 3 го; * кам пап пании ктора 1 года [
] ля НИИ / JJCd 0,0
эра да i>-
\ / //
/ /
// \й
/ ш. ^ /
с /1 /
'4 / / /
//
10и
133mXe/131mXe
10
Рис. 2. Смещение зоны реакторных выбросов в зону ядерных взрывов (сплошная линия — время кампании 3года, штриховая линия — время кампании 0,01 года)
парации изотопов Хе от предшественников и рассматривать все возможные гипотезы об источниках активности изотопов Хе.
Как было показано в работе [3], методы регуляризации являются одними из основных методов решения данной задачи. Измеренные значения активности изотопов записываются в виде системы линейных алгебраических уравнений (СЛАУ) АХ = Ь с плохо обусловленной матрицей АтА: а) для уравнений, содержащих активности изотопов и Кг, и Хе, матрица АтА имеет число обусловленности порядка 1018; б) для уравнений, учитывающих только активности изотопов Хе, число обусловленности уменьшается на три порядка. В настоящей работе решение СЛАУ проводилось методом регуляризации А.Н.Тихонова [8], а также исследовались другие методы: /р-регуляризация и векторная оптимизация, что особенно важно для решения задачи, когда измерены только изотопы Хе, где метод А.Н. Тихонова не всегда приводит к желаемому результату.
Кроме того, поскольку значения независимых выходов изотопов Кг и Хе имеют достаточно большие погрешности, то необходимо учитывать погрешности не только в правой части уравнений, но и в элементах матрицы А. Для учета погрешностей в элементах как матрицы Ь, так и матрицы А применяется конфлюэнтный анализ [6, 7].
В работе [3] предлагается способ определения времени сепарации изотопов Кг и Хе от их предшественников по цепочке радиоактивных превращений, без чего получить физически приемлемое решение невозможно.
Предлагаемый в работе [3] подход к решению задачи идентификации дает положительные результаты, но в ней не рассматриваются все возможные гипотезы формирования в атмосфере суммарной активности изотопов Кг и Хе, другими словами, не рассматриваются все возможные источники радиоактивных благородных газов.
Идентификация источников изотопов криптона и ксенона. Для надежной идентификации ядерных взрывов надо определить не только возможную принадлежность выброса к ядерным взрывам, реакторной активности, фоновой активности в атмосфере, но и определить возможный вклад в измеренную активность каждого изотопа некоторой постоянной составляющей (за счет возможной искусственной добавки активности в струе), т.е. набор гипотез значительно увеличивается. Задача будет считаться решенной, если можно на основе измеренных активностей изотопов благородных газов однозначно сказать, был ли ядерный взрыв или нет. Формально мы можем записать с учетом изотопов Кг и Хе не более 8-10 уравнений, определяющих активность измеренных изотопов, а для изотопов только Хе — не более 4 уравнений. Можно ввести некоторые модификации этих уравнений. Но число неизвестных будет намного больше, не менее числа
гипотез: три делящихся материала, две энергетические группы нейтронов, учет механизма сепарации изотопов Хе от предшественников, активация нейтронами продуктов деления, реакторная активность при различных временах кампании, вклад фоновой активности, вклад активности "заложенных" изотопов.
Система уравнений относительно неизвестных вкладов активностей изотопов благородных газов для различных видов деления записывается следующим образом (считаем, что время сепарации изотопов Кг и Хе от их предшественников по цепочкам радиоактивных превращений известно):
Ai
A2
X =
bi
b2
(1)
Строки матриц А1 и А2 ссоответствуют конкретному виду изотопов Хе и Кг, столбцы — активностям изотопов Хе и Кг для разных видов деления и для фона, элементы вектора Х показывают число делений данного вида, элементы векторов Ь1 и Ь2 — суммарные измеренные активности изотопов Хе и Кг соответственно.
Математические модели для расчета активности ¿-го изотопа для случая мгновенного деления и при функционировании реактора приведены [3].
Естественный природный фон можно приближенно рассчитать, однако далеко не всегда он обусловлен естественной природной активностью имеющихся в атмосфере изотопов Кг и Хе. Четких формул для расчета фоновой активности нет, поскольку ее наличие обусловлено совокупностью многих факторов, например, активностью "заложенных" изотопов. Но в математической модели для задачи идентификации ее необходимо учитывать.
Поскольку фоновая активность в общем случае является неизвестной, то в матрицах А1, А2 элементы, соответствующие фоновой активности изотопов Кг и Хе, примем соизмеримыми с активностями изотопов для разных видов деления (того же порядка).
Интервал, к которому принадлежит время сепарации, можно оценить на основе несложных графических построений. На рис. 3 представлена зависимость относительной активности изотопа 133тХе (в сравнении с активностью изотопа 135Хе) от времени. Относительная активность показана для двух видов деления (235И[1] и 239Ри[1\]) при мгновенном делении в зависимости от времени и в предположении, что изотопы не сепарируются от предшествующих им. Между линиями, соответствующими 235И[1] и 239Ри[1] должны проходить графики зависимости относительной активности от времени для всех остальных возможных видов мгновенного деления (235и[14], 238и[14],
Рис. 3. Изменение относительной активности изотопов 133mXe в случае деления 235ЩА (сплошная линия) и 239Pu[f] (штриховая линия) с учетом и без учета сепарации от предшествующих изотопов (точка 1 — отношение измеренных активностей изотопов в момент времени £ = 12 ч)
239Ри[14]), но чтобы не усложнять рис.3, они не показаны (символ в квадратных скобках указывает энергию нейтронов: / — энергия нейтронов спектра деления, 14 — нейтроны с энергией 14МэВ), т.е. показаны "верхняя" и "нижняя" линии, между которыми находятся графики относительной активности для других видов мгновенного деления — 235и[14], 238и[14], 239Ри[14].
Видно, что точка 1 (измеренная относительная активность на момент 12 ч) не принадлежит возможным значениям относительной активности изотопа без сепарации. Поэтому, чтобы оценить временной интервал ДЬ появления сепарации, рассчитаем относительную активность изотопа "в обратном времени", исходя из измеренных данных, в предположении, что отсутствует какое-либо влияние предшествующих изотопов на результаты измерения, т.е. активность изменяется по экспоненциальному закону (после сепарации изотопы распадаются по своим постоянным распада). Пересечение границ возможной относительной активности измеренного изотопа при мгновенном делении прямыми линиями показывает временной интервал ДЬ, когда произошло отделение изотопов от предшествующих им. Для рис. 3 он примерно равен одному часу.
Иногда качественный вывод о видах деления (делящегося материала и энергии нейтронов) можно дать уже исходя из подобных графиков. На рис.4 представлена зависимость относительной активности изотопов 85тКг и 88Кг (в сравнении с активностью изотопа 135Хе).
Рис. 4. Изменение относительной активности изотопов (сплошные линии)
и 88& (штриховые линии) с учетом и без учета сепарации от предшествующих изотопов (точки 1 и 2 — отношения измеренных активностей изотопов в момент времени £ = 6 ч)
Верхняя сплошная линия соответствует делению 235а нижняя — 239РиИ. Точки 1 и 2 — это отношения измеренных активностей в момент времени £ = 6 ч. Если аналогично рис. 3 рассчитать, начиная из точек 1 и 2, относительную активность изотопов "в обратном времени", то видно, что они пересекают соответствующие кривые (графики отношений без учета сепарации) в одно и то же время (почти в 3 ч). Отсюда следует, что образование изотопов благородных газов обусловлено в основном только одним видом деления. Если бы точки пересечения не совпали во времени, то это бы свидетельствовало о наличии нескольких видов делений.
При решении задачи идентификации источников изотопов благородных газов следует рассмотреть следующие причины образования радиоактивных изотопов в атмосфере:
1) реакторная активность и фоновая активность;
2) ядерный выброс (урановая бомба) и фоновая активность;
3) ядерный выброс (плутониевая бомба) и фоновая активность;
4) ядерный выброс (уран-плутониевая бомба) и фоновая активность;
5) ядерный выброс (урановая бомба), реакторная и фоновая активности;
6) ядерный выброс (плутониевая бомба), реакторная и фоновая активности;
7) ядерный выброс (уран-плутониевая бомба), реакторная и фоновая активности и др.
Сначала исследуем решение задачи по приведенным гипотезам различными математическими методами в предположении, что элементы матрицы А известны точно. Далее выберем те из них, которые дадут положительный результат, и применим для случая, когда элементы и матрицы А, и матрицы Ь заданы с известными погрешностями.
Поскольку число обусловленности матрицы Ат А имеет порядок 1018, воспользуемся методами регуляризации.
Методы регуляризации. Рассмотрим регуляризацию по А.Н. Тихонову и /р-регуляризацию [8, 9]. Для регуляризации Тихонова функционал имеет вид
7т (X) = || АХ - Ь||2 + Л ||Х||2; (2)
для /р-регуляризации
(X) = || АХ - Ь||2 + Л ||Х||р , 0 < р < 1. (3)
т
Здесь ||у| = ^^ |у^ |г; т — число компонент вектора у; Л — параметр
^=1
регуляризации (неотрицательное число), который подлежит дополнительному определению, причем четко формализованных процедур его вычисления в настоящее время нет и Л можно найти в результате калибровочных измерений.
Как правило, сглаживающий функционал ||х||2 (вторую норму) применяют, когда известно, что решение должно быть гладким. Если решение х будет иметь только часть ненулевых координат, а остальные равны нулю (или значительно меньше первых), то целесообразнее использовать сглаживающий функционал вида ||х||р, где 0 < р ^ 1.
Методы векторной оптимизации. Известно [8, 10], что методы регуляризации являются одним из способов сведения задачи векторной оптимизации к одной целевой функции с помощью параметра регуляризации Л. Здесь будем рассматривать исходную задачу как задачу векторной оптимизации и решать ее другими способами (метод е-ограничений, целевое программирование) без привлечения параметра регуляризации Л.
В методах регуляризации требуется найти минимум регуляризиру-ющего функционала 7 (X), который является суммой двух функционалов: функционала метода наименьших квадратов (||АХ — Ь||2) и сглаживающего функционала (||Х||2 или ||Х||р), объединенных параметром регуляризации Л, т.е. методы регуляризации можно рассматривать как частный случай метода взвешенных сумм.
Задачу векторной оптимизации запишем следующим образом:
J1 = II AX - b||2 ^ min, (4)
2 X
J2 = ||X||2 —^ min или J2 = ||X||« —^ min
2 X p X
при ограничениях
AX = b.
Ограничения в виде равенства AX = b можно было бы использовать в случае, когда элементы матрицы b известны без погрешностей. Поскольку на элементы вектора измерений b действуют помехи (мы будем рассматривать белый шум), то система ограничений-равенств в общем случае несовместна и равенства следует заменить на ограничения в виде двойных неравенств:
AX = b ^ b - 3Г < AX < b + 3Г,
где (Г = [ а1 а2 • • • а8 ]T — вектор средних квадратических отклонений (СКО) шума.
В целях упрощения записи функционалы J2 = ||X||2 и J2 = ||X||p будем записывать одним функционалом J2 = ||X||p, где P Е (0, 1]и[2].
Метод е-ограничений переводит все целевые функции, кроме одной, в ограничения. Применив этот метод к задаче (4), получим следующие виды задач математического программирования:
J1 = ||AX - b||2 ^ min при ||X||p < S, b - 3Г < AX < b + 3Г;
X (5)
J2 = ||X||p ^ min при ||AX - b|2 < e, b - 3Г < AX < b + 3Г.
X (6) В качестве оценки S можно принять максимальный физически возможный уровень активности либо максимальное значение, на которое рассчитаны датчики станций, регистрирующих радиоактивное излу-
[n~
чение; e можно принять равным 3* аг2, где N — число уравнений
у i=1
системы (1).
В целевом программировании существует две модели (два способа) решения задач векторной оптимизации: архимедова модель и модель с приоритетами. В архимедовой модели все целевые функции переводятся в ограничения и выполняется максимизация взвешенной суммы отклонений значений целевых функций от правой части ограничений. Для задачи (4) архимедова модель записывается следующим образом:
w1d1 + w2d2 ^ max, (7)
di, d2
где + = 1 при ограничениях
71 (Х) + й < е,
72 (Х) + ¿2 <
Ь — 3<х < АХ < Ь + 3<х.
В модели с приоритетами осуществляется последовательный перевод целевых функций в ограничения и максимизация отклонений значений целевых функций от ограничений. При этом найденное на данном шаге значение отклонения ^ используется как оптимальное отклонение на следующем г + 1 шаге.
Шаг 1.
тах при 71 (Х) + ^ е, Ь — 3сг ^ АХ ^ Ь + 3сг.
<1г
Шаг 2.
тах d2 при Ji (X) + dionx|dloni=dl = ^
J2 (X) + d2 < b - 3<г < AX < b + 3<r.
Линейное программирование. Для решения системы уравнений (1) можно также использовать линейное программирование (как частный случай при p = 1 и X ^ 0): минимизировать линейный функционал
F (X) = CTX ^ min (8)
при ограничениях
b - 3<г < AX < b + 3<г,
где C — вектор-столбец, состоящий из единиц.
Особенности идентификации источника по измеренной активности только изотопов ксенона. При идентификации источников изотопов Xe имеют дело с аналогичной (1) системой линейных уравнений. Но в этом случае она содержит только подматрицу A1 и вектор b1, т.е.
AiX = bi. (9)
При этом число элементов вектора X уменьшается до 11 — отсутствуют компоненты, соответствующие активностям изотопов Kr.
Результаты моделирования. В табл. 1 приведены результаты решения задачи идентификации ядерного взрыва для некоторых гипотез, полученные различными методами и при разных условиях (идентификация по изотопам Kr (83mKr, 85mKr, 85Kr, 87)Kr и Xe (131mXe, 133mXe, 133Xe, 135)Xe и идентификация только по изотопам Xe). Значения СКО решения получены методом статистических испытаний при действии
на компоненты вектора Ь белого шума с СКО, равным 1 % их номинального (точного) значения. Принималось, что момент сепарации равен 2 ч, момент измерений 24 ч, время кампании реактора 3 года.
Таблица 1
Гипотеза Элементы, Метод иден- Точное Найденное СКО
по изотопам тификации решение решение реше-
которых ния
проводилась
идентифика-
ция
Взрыв плу- Xe 1р -регуляриза- 100 60,5527 13,8580
тониевой ция 100 134,6635 17,8352
бомбы (А = 10-4, 0 2,6729 0,8593
Р = 0, 9) 0 0,0005 18,0917
0 -0,5955 2,6673
0 0,7365 2,4416
Реакторный Xe, Kr Линейное 100 100,0000 9,5230
выброс и программиро- 0 0 0
взрыв плу- вание 0 0 10,2339
тониевой 0 0.0000 16,1254
бомбы 0 0,0000 20,4340
100 100,0000 14,4489
100 100,0000 7,0458
0 0 0
0 0,0000 0
0 0,0000 0
0 0,0000 0
0 0 0
0 0 0
0 0 0
0 0 0
Реакторный Xe, Kr Регуляризация 100 108,1940 5,0176
выброс и А.Н. Тихонова 100 99,2889 13,6162
взрыв уран- (А = 10-5) 100 95,4071 7,1596
плутониевой 100 94,0900 6,4607
бомбы 100 101,5592 6,0731
Из всех методов, примененных для решения задачи (1), наилучшие результаты дали /р-регуляризация (при Л = 10-4, р = 0, 9) и линейное программирование, которые позволили получить верное решение для всех гипотез. При этом для точных исходных данных линейное программирование дает точное решение, в то время как результат решения методом /р-регуляризации несколько отличается от точного, но качественно верный.
Для решения системы (9) были применены те же подходы, что и для системы (1). Линейное программирование, позволяющее решить
задачу (1) при учете активностей изотопов Хе и Кг, при рассмотрении наиболее полной, 7-й и ряда других гипотез оказалось бесполезным — однозначного ответа, был взрыв или нет, дать нельзя. Решение удается получить только для более простых гипотез (1-й-6-й) методом /р-регуляризации (при Л = 10-4, р = 0, 9). Причина этого заключается в том, что, убирая из (1) уравнения с изотопами Кг, мы снижаем количество информации, которое может помочь в поиске решения.
2. Идентификация по результатам сейсмических измерений. Дополнительным критерием, позволяющим дать информацию о возможном проведении ядерных испытаний, является определение района испытаний. Ядерный взрыв (взрывы) можно попытаться провести скрытно, завуалировав его (на фоне землетрясения или другого взрыва). При этом, как правило, амплитуда сейсмических колебаний земной коры от ядерного взрыва в 10-15 раз меньше амплитуды от землетрясения. Поскольку сигнал сейсмических измерений колебаний земной коры, зарегистрированный при землетрясении, имеет широкополосный спектр, то обнаружение, а точнее пеленгация места проведения в это время ядерного испытания является более трудной задачей, нежели пеленгация источников гармонических сигналов.
Методы пеленгации источников сейсмической активности, в том числе в случаях попыток их сокрытия на фоне более мощного сейсмического события, описаны в [11]. Метод [11] позволяет решить задачу пеленгации в частотной области, откуда следует серьезный недостаток метода — для получения приемлемой точности решения в частотной области необходим большой объем исходной информации.
Учтем еще одну особенность задачи: регистрация сейсмических сигналов ведется непрерывно, но в некоторый момент времени в этой записи должен появиться сигнал от ядерного взрыва, визуально обнаружить который нельзя. Поэтому обработку сигналов необходимо вести непрерывно, разбив запись на временные фрагменты такой длительности, чтобы не пропустить запись ядерного взрыва.
Пусть имеется К источников сейсмической активности, различных по своей природе и имеющих разные мощности и формы сигналов. Отметим, что сигналы источников являются в общем случае негармоническими. Пеленгацию осуществляем посредством линейной системы сейсмических датчиков. За линию отсчета пеленгов примем линию, соединяющую датчики системы. Введем обозначения: в = [ в\ #2 ... #к ] — вектор пеленгов источников сейсмической активности; X = [ х\ х2 ... хк ]Т — вектор амплитуд (мощностей) излучаемых сигналов; Ь = [ Ь\ Ь2 ... Ьм ]Т — комплексная огибающая выходов датчиков системы, где М — число датчиков.
Ставится задача определения количества источников сейсмической активности, амплитуды (мощности) излучаемых ими сигналов, пеленгов источников.
Поскольку при распространении сигнала в земной коре на него накладывается помеха, а также имеют место ошибки измерений, обусловленные используемой аппаратурой, необходимо получить не только точечные оценки искомых параметров, но и их ковариационные матрицы или, по крайней мере, дисперсии. Помеху считаем аддитивной, имеющей нулевое математическое ожидание и известную дисперсию.
Основная идея решения задачи сейсмической пеленгации состоит в том, что посредством применения преобразования Фурье сигналы произвольной формы можно представить в виде суммы гармонических сигналов. Решение задачи пеленгации гармонических сигналов описано в работах [10, 12].
Все сигналы имеют одинаковые (близкие) спектры. Сначала предположим, что все источники сейсмической активности (их количество мы не знаем) излучают сигналы, имеющие одинаковые или близкие спектры. В данном случае можно применить следующий алгоритм решения задачи.
1. Проводим синхронную запись сигналов с выходов всех датчиков системы.
2. По выборке, полученной с базового датчика (крайний датчик), строим энергетический и фазовый спектры принятого сигнала. Нормируем энергетический спектр к единице.
3. Для энергетического спектра задаем шумовой порог Z — минимальную мощность полезных гармоник. Гармоники, по мощности меньшие порога, считаем шумовыми и при дальнейших вычислениях не учитываем.
4. Формируем матрицу А размера М х К, элементы которой вычисляются по формуле
N Г
Ат,к = ^2 ехр < у
9=1
2nf
2nfqt + (m - 1) dcos 6k + (fiq
где ад — относительная мощность д-й частотной составляющей, определенная по нормированному энергетическому спектру, ^ < ад ^ 1; ^ — частота д-й частотной составляющей, определенная по нормированному энергетическому спектру; — начальная фаза д-й частотной составляющей, определенная по фазовому спектру; N — число сигнальных гармоник (гармоник, превышающих по мощности шумовой порог); й — расстояние между соседними датчиками системы; V — скорость распространения сейсмических сигналов. Матрицу А можно
формировать с учетом нескольких временных отсчетов, полученных с датчиков системы [10].
5. Решаем СЛАУ АХ = Ь методами /р-регуляризации и векторной оптимизации, где Х — вектор распределения амплитуд (мощностей) источников сейсмической активности по пеленгам,.
Сигналы имеют произвольные спектры. Теперь рассмотрим общий случай, когда мы не знаем ни числа источников, ни спектров излучаемых ими сигналов. В данном случае можно применить следующий алгоритм.
1. Проводим синхронную запись сигналов с выходов всех датчиков системы.
2. По выборке, полученной с базового датчика, строим энергетический спектр принятого сигнала. Нормируем его к единице.
3. Для энергетического спектра задаем шумовой порог Z — минимальную мощность полезных гармоник. Гармоники, по мощности меньшие порога, считаем шумовыми и при дальнейших вычислениях не учитываем. Определяем число сигнальных гармоник.
4. Для каждой сигнальной гармоники решаем задачу одночастот-ной многосигнальной пеленгации. Формируем матрицу А размера М х К с элементами
Am,k = exp < j
2П/
2/t + (m - 1) d—cos 9k + (
где fq — частота q-й частотной составляющей, определенная по нормированному энергетическому спектру; (q — начальная фаза q-й частотной составляющей, определенная по фазовому спектру, полученному с базового датчика системы. Матрицу A можно формировать с учетом нескольких временных отсчетов, полученных с датчиков системы [10].
5. Формируем вектор b как виртуальную выборку временных комплексных отсчетов, полученную синхронно с выходов всех датчиков и представляющую собой комплексную амплитуду одночастотного сигнала частоты fq. Каждый элемент вектора b формируется по формуле aq exp {j [2n/qt + (q]}, где aq — относительная мощность q-й частотной составляющей, определенная по нормированному энергетическому спектру; (q — начальная фаза q-й частотной составляющей, определенная по фазовому спектру, полученному с соответствующего датчика системы (данная начальная фаза учитывает фазовый сдвиг, обусловленный геометрией системы датчиков). Решаем СЛАУ AX = b методами /р-регуляризации и векторной оптимизации, где X — вектор распределения амплитуд (мощностей) источников сейсмической активности по пеленгам.
Пример. Пусть имеют место два сейсмических события: подземный ядерный взрыв (азимут по отношению к системе датчиков 32,9°,
1,2 1,0
°'2 0 20 40 60 80 100 120 140 160 180
Рис. 5. Пеленгационная панорама
скорость распространения волны в земной коре 10,4 км/с) и землетрясение (азимут 104,4°, скорость волны 14,8 км/с). Предположим, что мощность взрыва приблизительно на два порядка меньше мощности землетрясения. Имитируется ситуация, когда взрыв хотят скрыть на фоне более мощного сейсмического события. Сначала рассмотрим ситуацию, когда на момент начала записи сигналов имеют место и взрыв, и землетрясение (рис. 5).
Для этого наложим сигналы от взрыва и от землетрясения друг на друга так, чтобы время взрыва совпадало с временем начала записи данных. Пеленг взрыва определен достаточно точно (32°). Пеленг землетрясения определен с большой ошибкой (148° вместо 104,4°). Это связано с тем, что математическая модель, рассмотренная выше, не адекватна реальным сейсмическим процессам. Неадекватность обусловлена тем, что в предложенной выше модели предполагается, что обе зарегистрированные системой датчиков сейсмические волны распространяются с одинаковой скоростью V. Таким образом, в данном случае был осуществлен расчет с подстановкой в модель скорости волны от взрыва. Поэтому его пеленг определен весьма точно. Для определения пеленга землетрясения воспользуемся тем, что скорость его сейсмической волны нам априори известна. Зная смещенный пеленг (148°) и реальную скорость распространения волны от землетрясения, учитывая, что косинус пеленга обратно пропорционален скорости распространения волны, проведем корректирующие вычисления и получим пеленг 104,8°. Таким образом, получили пеленги эпицентров обоих сейсмических событий, несмотря на то, что их сейсмические волны пришли на систему датчиков одновременно и мощность скрываемого события (в данном случае взрыва) на два порядка меньше мощности фонового события (в данном случае землетрясения).
Если в некоторый момент времени от начала записи сигналов системой датчиков имело место лишь фоновое сейсмическое событие, а спустя некоторое время, имел место взрыв, то алгоритм определит наличие только одного источника сейсмической активности — землетрясения, причем визуально по форме записанных с датчиков сигналов невозможно определить момент наступления скрываемого события. Это связано с тем, что, как правило, скрываемое событие имеет во много раз меньшую мощность, чем фоновое.
Таким образом, если разбить всю запись сигнала с системы датчиков на небольшие фрагменты и обработать каждый из них в отдельности, можно зафиксировать момент наступления скрываемого события.
Разработанный метод, позволяет осуществлять разделение и пеленгацию нескольких одновременно активных эпицентров сейсмических событий. В общем случае метод позволяет осуществлять многосигнальную пеленгацию источников негармонических сигналов.
Учет неопределенностей элементов матрицы А вместе с неопределенной правой частью Ь (конфлюэнтный анализ). В реальной ситуации помимо того, что измеренные значения Ь известны с ошибкой, элементы матрицы А также заданы неточно. Это обусловлено тем, что, например, в задаче идентификации источников изотопов благородных газов независимые выходы продуктов деления, которые используются при расчете коэффициентов матрицы А, имеют большую погрешность. Для идентификации по результатам сейсмических измерений одной из причин наличия погрешности в элементах матрицы А является неточность определения расстояния между сейсмическими станциями и скорости распространения колебаний в земной коре. Поэтому в алгоритмах получения решения Х необходимо предусмотреть процедуры учета всех погрешностей.
Рассмотрим систему уравнений
АХ = Ь (10)
в общем случае (для идентификации по изотопам и сейсмике).
Система (10) — линейная относительно неизвестных х^, поэтому рассмотрим возможность ее решения в предположении, что элементы матрицы А системы и правой части Ь заданы с погрешностями, распределенными по нормальному закону с математическими ожиданиями, равными нулю, и дисперсиями, равными соответственно а2 (а*™6) и а2 (Ь*гае):
агз = а?" ± , - N (0, а2 )) , Ь = Г ± й, й - N (0, а2 (Се)) ,
где b
^rue i '
7true
Hj
— истинные значения активностей; N обозначает нор-
мальное распределение.
Для учета погрешностей как в правой части, так и в элементах матрицы системы используем определение ортогональной регрессии [13]:
( ( т \ 2 \
i=1
bi-E
aij xj
о2 (birue)
+
m
E
j=1
Hj
-atrue)2
- ij
о2
airUe)
V
m
+a> :ixjip. (ii)
j=i
/
В функционале (11) принимается, что ошибки — статистически независимые величины. Наряду с неизвестным вектором х^- в функционале (11) неизвестны также истинные значения вычисляемых активностей а*™6. Единственность и состоятельность оценок, получаемых методами конфлюэнтного анализа, доказаны в работе [7].
Структурная схема алгоритма приведена на рис. 6.
На первом шаге алгоритма положим а*™6 = а^ и, решив задачу (11), найдем решение х^.
Для получения оценок истинных значений а у16 (в случае идентификации по изотопам — при заданном значении ¿0, по сейсмике — при заданных d и V) на каждом шаге вычисления оценок х^- используется условие
арк = 0,
daje
что приводит к решению систем линеиных уравнении следующего вида:
true j true
XpO
E
=i о2 (6Г)
+
a.
ip
о2 la.
true
''ip J
+
о2 №) о2 (ЬГ )■
Полученные новые оценки значений а*'"6 должны удовлетворять естественному условию (принадлежать области неопределенности измеренных значений а^):
к, - а!™6| < 3а (а*™6) .
ij
ij
ij true
Если это условие не выполняется, то а*™6, которые не удовлетворяют этому неравенству, следует заменить на значения ближайших граничных точек. Из-за этого может происходить увеличение значений функционала при новых точных значениях переменных по сравнению с предыдущим шагом итерационного процесса, что приводит к снижению скорости сходимости итерационного процесса и возникновению
true
true
ЧГ=аГ
Рис. 6. Схема алгоритма решения плохо обусловленной СЛАУ с помощью кон-флюэнтного анализа
колебаний. Чтобы значения функционала не увеличивались после пересчета оценок а*™6, те наборы оценок а*™6, для которых произошло увеличение слагаемых функционала по сравнению с их значениями на предыдущей итерации, следует заменить на соответствующие значения для предыдущего шага.
После определения оценок истинных значений анаходим очередное приближение к решению х^, используя методы регуляризации.
Критерием останова алгоритма является несущественное различие значений функционала и компонент вектора Xj на соседних итерациях, т.е. выполнение неравенств
(xj)l - (xj У 1
(xj)
< Yi;
Fk ((Xj)l-1) - FK ((Xj)l)
Fk ((Xj)l)
< Y2,
где (х^) - очередное приближение к решению на 1-й итерации; 0 <71 < 1, 0 <72 < 1 — некоторые числа.
Для задачи идентификации ядерных взрывов по результатам измерения активности благородных газов существует еще один этап — определение времени сепарации изотопов Хе и Кг. Минимум функционала конфлюэнтного анализа по переменной ¿0 может быть найден каким-либо из методов одномерной минимизации. Значение ¿Т™, при котором функционал достигает минимума, будет искомым значением времени сепарации изотопов от предшествующих им.
Общая схема алгоритма, позволяющего найти оценки а0, х^, а*™6, приведена на рис. 7.
После того как получены точечные оценки решения, необходимо найти интервальные оценки решения (или как минимум СКО). Согласно работе [14] ковариационная матрица оценок определяется
соотношением
. д 2F
Dn = -
( ÖXiÖXj
1
(12)
где F — соответствующий функционал, определяющий оценку реше-
ния х^.
Таким образом, СКО оценок имеют вид
&гг = \fDii.
Аналитические выражения для производных функционала (11) довольно сложны, поэтому для вычисления ковариационной матрицы
Рис. 7. Схема алгоритма получения оценок to и решения x
(12) используем функционал ^, также учитывающий погрешности как в правой части, так и в матрице системы [13]:
т \ 2 Ь _ а3 ) т = £ —-£-'— + А £ |р,
г=1 а2 (6г) + £ а2 (а„) ж? з=1
3 = 1
полученный из функционала путем переноса ошибок из левой части системы (10) в правую.
Этот функционал эквивалентен функционалу (11). Для других методов (например, методов векторной оптимизации) после перехода к задаче математического программирования ковариационная матрица может быть получена с помощью функции Лагранжа.
Выводы. Рассмотрены два подхода для решения задачи идентификации ядерного взрыва (и соответствующие им математические модели): идентификация по измеренным в атмосфере активностям изотопов Кг и Хе и идентификация по результатам сейсмических измерений. Показано, что каждая задача требует своего способа решения (методы регуляризации А.Н. Тихонова, /р-регуляризации, векторной оптимизации). Получено, что идентификация источников радиоактивности только по изотопам Хе возможна для самых простых гипотез (реакторный выброс, ядерный взрыв урановой или плутониевой бомбы при наличии слабого фона) в силу плохой обусловленности матрицы системы. Идентификацию источников радиоактивных газов для всех возможных гипотез можно осуществить только при учете активностей как изотопов Хе, так и изотопов Кг. Необходимо обратить особое внимание на регистрацию изотопов Кг. Для этого достаточно в разрабатываемых станциях регистрации изотопов Хе вывести гамма-спектр газообразных продуктов деления, но не отсекать гамма-излучение изотопов Кг, как это имеет место в настоящее время.
Дополнительную информацию для принятия гипотезы о наличии ядерного взрыва предоставляет идентификация по результатам сейсмических измерений, что делает решение более достоверным.
Предлагаемые методы позволяют выявить "заложенные" изотопы. В процессе идентификации ядерного взрыва определяются точечные и интервальные оценки найденных характеристик.
СПИСОК ЛИТЕРАТУРЫ
1. Г р е ш и л о в А. А., К о л о б а ш к и н В. М., Д е м е н т ь е в С. И. Продукты мгновенного деления и235, и2зв, Ри2зэ в интервале 0-1 ч.: Справочник. -Атомиздат, 1969. - 104 с.
2. Грешилов А. А., Колобашкин В. М. Способ определения концентрации изотопов инертных газов в смеси продуктов деления. АС № 366771 с приоритетом от 01.11.1967 г.
3. Грешилов А. А. Тетюхин A. A. Алгоритм идентификации источников радиоактивных благородных газов // Вестник МГТУ им. Н.Э. Баумана. Сер. "Естественные науки". 2005. - № 6. - С.28-37.
4. Martin B. Kalinowski, Christoph Pistner // Journal of Environmental Radioactivity. - V. 88 (2006). - P. 215-235.
5. Гусев Н. Г., Рубцов П. М., Коваленко В. В., Колобашкин В. М. Радиационные характеристики продуктов деления: Справочник. - М.: Атомиздат, 1974. - 224 с.
6. Г р е ш и л о в А. А. Математические методы принятия решений: Учеб. пособие для вузов. - М.: Изд-во МГТУ им. Н.Э. Баумана, 2006. - 584 с.
7. Г р е ш и л о в А. А. Анализ и синтез стохастических систем. Параметрические модели и конфлюентный анализ. - М.: Радио и связь, 1990. - 320 с.
8. Тихонов А. Н., АрсенинВ. Я. Методы решения некорректных задач. -М.: Наука, 1979,- 142 с.
9. C e t i n M. and Karl W. C. Feature-enhanced synthetic aperture radar image formation based on nonquadratic regularization // IEEE Trans. Image Processing. -2001. - Vol. 10, no. 4. - P. 623-631.
10. Грешилов А. А., Лебедев А. Л., Плохута П. А. Многосигнальная пеленгация источников радиоизлучения на одной частоте как некорректная задача // Успехи современной радиоэлектроники. - 2008. - № 3. - C. 30-46.
11. Кушнир А. Ф. Оценивание вектора кажущейся медленности плоской волны по данным трехкомпонентной сейсмической группы: статистическая задача с мешающими параметрами // Теоретические проблемы в геофизике: Сб. науч. труд. - М.: Наука, 1997. - 235 с. (Вычислительная сейсмология. Вып. 29). -C. 197-214.
12. Г р е ш и л о в А. А., П л о х у т а П. А. Многосигнальная пеленгация источников радиоизлучения на одной частоте / Вопросы защиты информации: Науч.-практ. журн. / ФГУП "ВИМИ". - 2008. - Вып. 1(80). - C. 61-67.
13. Г р е ш и л о в А. А., С т а к у н В. А., С т а к у н А. А. Математические методы построения прогнозов. - М.: Радио и связь, 1997, - 112 с.
14. Грешилов А. А., С т а к у н В. А., С т а к у н А. А. Статистические методы принятия решений с элементами конфлюентного анализа. - М.: Радио и связь, 1998,- 112 с.
Статья поступила в редакцию 23.06.2008
Анатолий Антонович Грешилов родился в 1939 г., окончил в 1964 г. Московский инженерно-физический институт. Д-р техн. наук, профессор кафедры "Вычислительная математика и математическая физика" МГТУ им. Н.Э. Баумана. Автор более 150 научных работ, в том числе 24 монографий, 20 авторских свидетельств и патентов в области разработки математических методов строгого учета неопределенности исходной информации в задачах математической физики, распознавания образов, прогнозирования и других технических приложений.
А.А. Greshilov (b. 1939) graduated from the Moscow Institute for Engineering and Physics in 1964. D. Sc. (Eng.) professor of "Computational Mathematics and Mathematical Physics" department of the Bauman Moscow State Technical University. Author of more than 150 publications including 24 monographs, 20 author's certificates and patents in the field of development of mathematical methods for strict account of source data uncertainty in problems of mathematical physics, image identification, forecast and other technical applications.
Алексей Леонидович Лебедев родился в 1983 г., окончил в 2007 г. МГТУ им. Н.Э. Баумана. Аспирант кафедры "Высшая математика" МГТУ им. Н.Э. Баумана. Специализиуется в области методов решения некорректных задач.
A.L. Lebedev (b. 1983) graduated from the Bauman Moscow State Technical University in 2007. Post-graduate of "Higher Mathematics" department of the Bauman Moscow State Technical University. Specializes in the field of methods of solving ill-posed problems.
Павел Анатольевич Плохута родился в 1983 г., окончил МГТУ им. Н.Э. Баумана в 2006 г. Аспирант кафедры "Вычислительная математика и математическая физика" МГТУ им. Н.Э. Баумана. Автор четырех научных работ в области решения некорректных задач и радиопеленгации.
P.A. Plokhuta (b. 1983) graduated from the Bauman Moscow State Technical University in 2006. Post-graduate of "Computational Mathematics and Mathematical Physics" department of the Bauman Moscow State Technical University. Author of 4 publication in the field of solving ill-posed problems in radio-direction finding.
