Алгоритм расчета плотности потока радона (222Rn) с поверхности земли Текст научной статьи по специальности «Математика»

ВЕСТНИК ТОМСКОГО ГОСУДАРСТВЕННОГО УНИВЕРСИТЕТА

2008 Математика и механика № 3(4)

УДК 550.338

Р.И. Паровик, П.П. Фирстов

АЛГОРИТМ РАСЧЕТА ПЛОТНОСТИ ПОТОКА РАДОНА (22^п)

С ПОВЕРХНОСТИ ЗЕМЛИ

Предложен алгоритм расчета плотности потока радона (ППР), на основе которого разработана программа («РЭКСЭМ») его расчета и визуализации одновременно с исходными данными. Программа прошла тестирования на данных, полученных сетью станций мониторинга подпочвенного радона в районе Петропавловск - Камчатского геодинамического полигона.

Ключевые слова: радон, плотность потока радона, перенос радона.

В последние десятилетия уделяется повышенное внимание к поиску предвестников землетрясений в сейсмоэманационных эффектах. В работах [1 - 4] показана перспективность таких исследований и предложен ряд новых подходов, в которых рассматриваются методы исследования массопереноса почвенного радона, с целью поиска предвестников сильных землетрясений, причем расчеты проводятся с помощью диффузионно-конвективной модели [5]. На Петропавловск-Камчатском геодинамическом полигоне регистрация Rn в каждом пункте осуществляется в двух точках на различных глубинах в рыхлых отложениях, которые рассматриваются как относительно однородная и пористая среда [6]. Как показано в работе [4], для такой среды в динамике Яп будет слабо выражена реакция на изменение напряженно-деформированного состояния геосреды. Это связано с тем, что изменение скорости конвекции не может привести к росту объемной активности (ОА) Яп выше ее фонового значения [7].

Авторы работы [7] предлагают для повышения эффективности выделения предвестниковых аномалий использовать параметр, связанный с ростом скорости конвекции - плотность потока радона (ППР) с поверхности земли. Увеличение скорости конвекции приводит к увеличению конвективного потока и как следствие к увеличению ППР на земной поверхности. Расчеты показали, что увеличение скорости конвекции в два раза приводят к росту ППР в 4 раза, в тоже время рост концентрации Яп увеличивается всего в два раза, причем в неоднородной среде ППР ещё более чувствительна к вариациям скорости конвекции. В работе [8] была реализована эта идея, где расчет ППР с поверхности земли проводился с помощью эмпирической формулы [9].

Естественно, что при поисках предвестниковых аномалий необходимо учитывать влияние вариаций метеорологических величин, в частности атмосферного давления. Исходные данные, полученные на сети станций, подвергались барокомпенсации по методике А.А. Любушина.

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

значения в зависимости от превышения над уровнем фона. Данная методика прошла тестирование на данных сети станций мониторинга подпочвенного радона на Камчатке [10].

1. Вычисление плотности потока на земной поверхности как решение обратной задачи массопереноса радона

1.1. Постановка задачи

Задача нахождения ППР на земной поверхности по измеренным значениям концентраций радона Ыех на различных глубинах с краевыми условиями имеет следующий вид [5] :

дЫ (г, < ) = в _ ^«(£,0 _ш г> , ) + е>

дї е дг2 дг

N(2,0) = N^, г = 0,

дN (г, г) (1)

= ч(г),

дг

.Ы (, 1 ] ) = Ыех ■

Здесь N = Q/n^, Бк/м3 - фоновая концентрация радона на заданной глубине для данного района мониторинга; 2, Бк/(м3-с) - скорость образования (эманирования) радона в грунте; п- пористость грунта; ^, с-1 - постоянная полураспада радона; Д м2/с - коэффициент диффузии радона в грунте; Д=Д/п, м2/с - эффективный коэффициент диффузии; V, м/с - скорость конвективного потока радона; qc, Бк/(м2-с) - постоянная единичная ППР на земной поверхности; М^,/), Бк/м3 - объемная концентрация радона. В данной модели предполагается, что процесс массо-переноса радона на земную поверхность в некоторый момент времени является установившимся. Поэтому параметры модели являются константами и также предполагаются известными.

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

1.2. Алгоритм решения обратной задачи массопереноса радона

Рассмотрим алгоритм решения задачи по определению ППР, основываясь на приближенном соответствии между измеренной концентрацией Яп для нескольких разноглубинных датчиков и рассчитанной концентрации по диффузионноконвективной модели [5,12].

Алгоритм решения можно разделить на следующие этапы:

1) аппроксимация искомой ППР постоянной единичной ППР и решение соответствующей прямой задачи по определению концентрации Яп в точках наблюдения;

2) представление решения прямой задачи массопереноса Яп как суперпозиции его элементарных блоков с помощью интегральной формулы Дюамеля;

3) минимизация невязки между концентрацией Яп, полученной экспериментальным путем и рассчитанной по модели методом наименьших квадратов.

1.3. Аппроксимация ППР ступенчатой единичной функцией

Для решения прямой задачи массопереноса радона аппроксимируем ППР на земной поверхности постоянной единичной функцией:

На поверхности ППР аппроксимируется ступенчатой единичной плотностью. На интервалах времени т1/2, т3/2,...., тМ-1/2 используются значения плотности потока qbq2,..., qM в моменты времени ім: от 0 до т1, от т1 до т2, от тМ-1 до тм.

1.4. Решение прямой задачи массопереноса радона

Прямая задача массопереноса радона на земную поверхность представляет собой дифференциальное уравнение в частных производных с начальными и граничными условиями [І2]:

Это решение описывает процесс изменения концентрации Яп в рыхлых отложениях под воздействием постоянной единичной плотности потока qc в моменты времени г. Дифференцируя решение (3) по qc, получаем выражение для коэффициента чувствительности ф(г,г) = 5ЛГ(г,г)^с = - да2ДеН/(^п), причем если коэффициенты чувствительности малы или коррелированны, то задача оценивания становится трудной и чувствительной к погрешностям измерений.

Согласно аппроксимации ППР в виде единичной функции и теореме Дюамеля о представлении решения в виде суперпозиции элементарных блоков его решений по соответствующим интервалам времени имеем

Д(г,/М) = Д+д^ф^м-То) - ф^м - Ті)]+д2[ф(2,гм - Ті) - ф(г,?м - т2)]+....+

+ дм[ф(г,гм - Тм-і) - ф(2,?м - Тм)],

Здесь ф(г, гм) - элементарное решение в 1М момент времени. Выражение (4) можно упростить, записав его в следующей дискретной форме:

^ N(z,0) = No, t = О,

(2)

-De ndN ( ’1) + vnN (о, t ) = qc.

д z=0

Решение задачи (2) можно получить интегральным преобразованием:

N(z,t) = (N - qc) m2DeH/(^n),

H = H1 + H2 - H3, H1 = (m + n)exp((m - n)z) erfc(z/(4Det)1/2 - n(Dei)1/2), H2 = (m - n) exp((m + n)z) erfc(z/(4Det)1/2 + n(Dei)1/2),

H3 = 2exp(-zm - ^i) erfc(z/(4Det)1/2 - m(Det)1/2), m = v/(2De), n = (v2 + 4De^)1/2/(2De).

(З)

1.5. Интегральная формула Дюамеля

ф(^ ім - Тм) = фМ)], M = 0, 1, 2,..., i.

(4)

м

N (г, 1м ) = ^ + Е qnAy (г, 1м-„ ), Аф (г, 1м-„ ) = Ф (г, ^-п-и )-ф (г, 1м-„ ), (5)

П= 1

которое называется интегральной формулой Дюамеля. Опуская индексы его можно переписать в следующей простой форме:

м

КМ = + Х ЧпЛФм-п ■ (6)

п= 1

При решении обратных задач массопереноса радона на земную поверхность, выражение (6) является важным, так как устанавливает удобную связь между концентрацией Rn и ППР на поверхности. Так как ППР изменяется по времени, то выражение (6) дает численный результат.

1.6. Метод наименьших квадратов для обратной задачи массопереноса радона

Если у нас имеется J датчиков измерения концентрации Rn, расположенных на различных глубинах, то тогда для каждой точки наблюдения получаем:

К\Ы = ^1м|?м =0 ,

^гш = ^2м|?м =0 +^21^М, (7)

= -^уш|?м =0 +фу 1 Чш ■

Между измеренной и рассчитанной концентрациями Яп довольно сложно установить точное согласование, так как единственную составляющую ППР трудно выбрать чтобы выполнялось равенство Щм = -^е1. Поэтому согласование должно выполняться в некотором усредненном приближении. Оценим дм с учетом (7) методом наименьших квадратов:

J . .2

5 =

]=

Минимизируем функционал невязки (8) по составляющей плотности потока дм, что сводится к нахождению нулей производных по параметру дм:

J г М-1 л

X I Кех - X Я л Аф ум - - ^ 1ф л

Ям = —---------^----------------— . (9)

2

Е ф£1

к=1

Замечание. Предлагаемый алгоритм решения обратной задачи по нахождению ППР имеет один недостаток - чувствительность к погрешностям измерения при малых шагах по времени (секунды) [12]. Однако регистрация экспериментальных данных сети мониторинга радона осуществляется с шагом дискретизации 6,0 ч-1, поэтому при реализации данного алгоритма чувствительность к погрешностям измерения является не высокой.

2. Экспериментальные данные и результаты расчетов

В пункте «Институт» (ИНС), располагающемся рядом со зданием Института вулканологии и сейсмологии ДВО РАН, осуществляется регистрация объемной активности (ОА) Яп в подпочвенном воздухе на глубине один метр от пола под-

земного бункера глубиной 2,5 м (зона аэрации) и с поверхности пола. Над бункером на поверхности земли установлен металлический контейнер, причем подземный бункер через контейнер и систему труб вентилируется за счет естественной конвекции воздуха. В наземном контейнере и подземном бункере установлены датчики температуры и атмосферного давления. Регистрация всех параметров осуществляется с помощью двух измерительных приборов ALMEMO 2590-9 с частотой дискретизации 6,0 ч-1 [8].

Как показало сравнение метеорологических величин в подземном бункере и наземном контейнере, в подземном бункере в летнее время, за счет больших вариации температуры в наземном контейнере, четко прослеживается суточные колебания атмосферного давления, которые хорошо видны на рис. 1, а, где показана динамика ОА Rn за период 20 июля - 29 августа 2006 г. в зоне аэрации и с поверхности. На рис. 1, а видно, что с 13 по 22 августа наблюдается увеличение амплитуды объемной активности Яп на обоих каналах, причем в зоне аэрации превышение составляет 115%, а на поверхности 130% от фона. На рис. 1, б приведена расчетная кривая ППР, вычисленная по описанной выше методике с помощью программы «РЭКСЭМ» со следующими значениями параметров: А1 = 0,2 м, А1 = 1,2 м, Б = 5,5-10-7 м2/с, V = 5• 10-5 м/с, X = 2,1-Ю-6 с-1. В расчетных кривых ППР превышение амплитуды над фоном для аномального периода составляет более 160%, что подтверждает более высокую чувствительность ППР по сравнению с динамическими параметрами ОА Яп.

н

О

О

я

со

я

5

3

я

*0

ю

О

ш 'у'

о ^

к ^

* Ш

о

я

н

о

ч

С

Рис. 1. Динамика объемной активности радона (а) и вычисленная плотность потока радона (б) за период 20 июля - 29 августа 2006 г. 1 - зона аэрации, 2 - поверхность пола подземного бункера

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

ЛИТЕРАТУРА

1. Рудаков В.П. Сейсмоэманационные эффекты геологических структур // Проблемы геофизики XXI в. Кн. 2. М.: Наука, 2003. С. 95 - 113.

2. Dubinchuk V.T. Radon as a precursor of earthquakes // Isotopic geochemical precursors of earthquakes and volcanic eruption. Vienna, 1991. P. 6 - 22.

3. Virk H.S., Baljinder Singh. Radon anomalies in soil-gas and groundwater as earthquake precursor phenomena // Tectonophysics. 1993. V. 227. P. 215 - 224.

4. Nasaroff W.W., Nero A.V. (Eds). Soil as a source of indoor radon: generation, migration and entry. Radon and its decay products in indoor air // A Wiley-Interscience publications, 1988. P. 57-112.

5. Новиков Г.Ф., Капков Ю.Н. Радиоактивные методы разведки. М.: Недра, 1965. 750 с.

6. Фирстов П.П., Рудаков В.П. Результаты регистрации подпочвенного радона в 1997 -2000 гг. на Петропавловск-Камчатском геодинамическом полигоне // Вулканология и сейсмология. 2003. № 1. С. 26 - 41.

7. Яковлева В.С, Каратаев В.Д. Плотность потока радона с поверхности земли как возможный индикатор изменений напряженно-деформированного состояния геологической среды // Вулканология и сейсмология. 2007. № 1. С. 74 - 77.

8. Фирстов П.П., Широков В.А., Руленко О.П., Яковлева В.С, Исаев А.В., Малышева О.П. О связи динамики подпочвенного радона (222Rn) и водорода с сейсмической активностью Камчатки в июле-августе 2004 г. // Вулканология и сейсмология. 2006. № 5. C. 49 - 59.

9. Рыжакова Н.К., Яковлева В.С. Патент РФ № 2212688 от 20.09.2003. Способ определения плотности потока с поверхности земли.

10. Фирстов П.П., Пономарев ЕА., Чернева Н.В., Паровик Р.И. Исследование кинематических и динамических параметров эманации подпочвенного радона в период активации сейсмичности Камчатки в августе 2006 г. // Солнечно-земные связи и предвестники землетрясений: Сб. науч. докл. IV Междунар. науч. конф., с. Паратунка, Камчатская обл., 2007 г. С. 464 - 469.

11. Паровик Р.И. Программа обработки геофизических данных «РЭКСЭМ» // Телеграф отрасли фонда алгоритмов и программ. Инновации в науке и образовании. М.: Издательский дом Святогор, 2007. № 5. С. 19.

12. Бек Дж., Блакуэлл Б., Сент-Клэр Ч. Некорректные обратные задачи теплопроводности. М.: Мир, 1989. 312 с.

Статья принята в печать 08.09.2008г.