CC BY
81
УДК 004.9:504:519.6
И.В. КОВАЛЕЦ, С.А. ЮЩЕНКО, М.И. ЖЕЛЕЗНЯК
АССИМИЛЯЦИЯ ДАННЫХ В ТРЕХМЕРНОЙ МОДЕЛИ РАСПРОСТРАНЕНИЯ РАДИОНУКЛИДОВ В ВОДОЕМАХ ТРИТОКС СИСТЕМЫ РОДОС_______________________________________________
Abstract: The data assimilation procedure is developed for the 3D model of the radionuclide dispersion in water systems THREETOX included into the Hydrological Disperision Module of the EC decision support system for offsite nuclear emergency management RODOS. The method is verified on the basis of the monitoring data on radioactive contamination of the Black Sea, and also on the simulated data of the radionuclide dispersion in the Kakhovka Resrvoir for the scenario of hypothetical accidental atmospheric fallout from the Zaporizhke NPP.
Key words: data assimilation, radionuclide dispersion, three-dimensional model.
Анотація: Розроблений алгоритм асиміляції даних вимірювань концентрацій радіонуклідів у тривимірній моделі розповсюдження радіонуклідів у водних системах ТРИТОКС, яка включена в гідрологічний модуль європейської системи підтримки прийняття рішень при ядерних аваріях РОДОС. Метод верифікований на основі даних вимірювань у Чорному морі, а також на основі результатів моделювання розповсюдження радіонуклідів у Каховському водосховищі для сценарію атмосферного випадання від гіпотетичного аварійного викиду із Запорізької АЕС.
Ключові слова: асиміляція даних, розповсюдження радіонуклідів, тривимірна модель.
Аннотация: Разработан алгоритм ассимиляции данных измерений концентраций радионуклидов в трехмерной модели распространения радионуклидов в водных системах ТРИТОКС, входящей в гидрологический модуль европейской системы поддержки принятия решений при ядерных авариях РОДОС. Метод верифицирован на основании данных измерений в Черном море, а также на основании результатов моделирования распространения радионуклидов в Каховском водохранилище для сценария атмосферных выпадений от гипотетического аварийного выброса из Запорожской АЭС.
Ключевые слова: ассимиляция данных, распространение радионуклидов, трехмерная модель.
1. Введение
Трехмерная модель ТРИТОКС [1] включена в гидрологический модуль европейской системы поддержки принятия решений по внеобъектному реагированию на радиационные аварии РОДОС [2] для описания миграции радионуклидов в глубоких водоемах с выраженной вертикальной изменчивостью гидродинамической структуры потока и полей загрязнений [3]. ТРИТОКС может начинать моделирование на основе сценариев поступления радионуклидов в водоем из притоков, рассчитанных одномерной речной молелью РИВТОКС [3] и/или на основе полей атмосферных выпадений на поверхность водоема, рассчитанных моделями атмосферного модуля системы РОДОС на основе сценариев аварийных выбросов из реактора АЭС. С течением некоторого времени после выброса появляются данные измерений концентраций радионуклидов, которые могут использоваться для уточнения прогнозов распространения радионуклидов в водоеме в последующие моменты времени. Проблема эффективного использования данных измерений для повышения точности прогностических моделей известна как проблема ассимиляции данных или усвоения данных [4]. Опыт пилотного использования системы РОДОС показал необходимость развития методов ассимиляции данных в прогностических моделях, используемых в системе РОДОС [5].
Методы ассимиляции данных, применяемые в океанографии и метеорологии, достаточно полно описаны в [6]. Эти методы могут быть условно разделены на два класса: статистические и вариационные методы. Общим статистическим методом является Калмановская фильтрация [7], для которой совместно с системой прогностических уравнений основных расчетных величин (скорости, концентрации и других) должны решаться прогностические уравнения для компонент
ковариационных матриц. В вариационных методах решается задача оптимального управления с целью уточнения начальных и граничных условий, а также параметров модели [8].
В рамках проекта РДОЫЕМ ”Пятой рамочной программы ЕС“ были разработаны методы ассимиляции данных для прогностических моделей системы РОДОС, включая модели гидрологического модуля этой системы [4]. В связи с тем, что система должна работать в реальном времени, вычислительная эффективность является одним из основных требований, что существенно ограничивает выбор из существующего ряда методов ассимиляции данных. Это делает упомянутые выше подходы мало приемлемыми для трехмерной модели.
В данной работе представлен статистический метод ассимиляции данных, близкий к так называемой оптимальной интерполяции. В этом методе, в отличие от алгоритма Калмановской фильтрации, компоненты ковариационных матриц предполагаются заранее заданными. Представлены результаты верификации этого метода на основе моделирования распространения радионуклидов в Черном море и Каховском водохранилище.
2. Описание метода
ТРИТОКС является современной трехмерной моделью гидротермодинамики поверхностных вод и распространения радионуклидов в эстуариях, озерах и морях [1]. ТРИТОКС включает в себя подмодели гидродинамики, переноса взвешенных наносов и переноса радионуклидов. Гидротермодинамика описывается в рамках трехмерной нестационарной негидростатической модели. Ее прогностическими переменными являются три компоненты скорости, температура, соленость и возвышение свободной поверхности. Перенос наносов описывается адвективнодиффузионным уравнением с членом типа «источник-сток», описывающим взмучивание и осаждение взвесей. Изменение толщины верхнего слоя донных отложений описываются уравнением деформации дна. Уравнения переноса радионуклидов моделируют распространение радионуклидов в трех фазах: в растворе, на взвешенных наносах и в донных отложениях. Обмен между этими фазами представлен в уравнениях членами, описывающими процессы адсорбции-десорбции и взмучивания-осаждения. Уравнения модели ТРИТОКС интегрируются численно с помощью конечно-разностных методов.
Для ассимиляции данных в ТРИТОКС после тестирования различных методов был избран метод “итераций к оптимальному решению” (ИОР). Этот метод был впервые предложен в [6] и применялся в атмосферных моделях системы РОДОС [9]. Как было показано в [6] , метод ИОР
близок к известной “оптимальной интерполяции” [6]. В алгоритме ИОР поле концентрации С^ (X,/)
обновляется для времени t = tl в узлах расчетной области . = 1,N. в итерационном цикле:
б“ (л+1)( X, ,^) = б? (п)( X, ,^) +
(1)
+(1 + )-1[СГ (0)( X, Л) - б:(п)Х Л)];
(2)
где Х{ =(хи,х1{) - радиус-вектор точки .; п - количество итераций; б'1 (0)(X,tl)- поле первого
приближения концентрации перед началом итерационной процедуры; б'МоЬз(Xr,to) - измеренная концентрация радионуклидов в точках г = 1,..,Ы в момент времени t = tO, ближайший к моменту времени t = tl; сЦс20 - отношение среднеквадратического отклонения поля первого приближения к среднеквадратической ошибке наблюдений.
Весовая функция w1(Xk -X,) • - tO) = ,Stol) является корреляционной функцией
отклонений поля первого приближения от истинного поля, так что -Xi) = 1 при Xk = Xi, и
-Xi) ® 0 при ав Xк -X\[ . Аналогичные соотношения верны для -tO). В данной
работе используются соотношения
м1(X - X ) = еХР^^Т-Т),гГк = >/ (Х1к - ХИ )2 + (Х2к - Х2, )2
2Коог
2 , (3)
^ - Ъ ) = еХР(-2^4^ 3{1о = - {О
2Я1оог
где Яоог, Я1оог- пространственный и временной радиусы влияния.
Преимущества метода ИОР следующие: (1) метод практически не требует дополнительного к времени расчета прогноза распространения концентрации вычислительного времени, что чрезвычайно важно для систем реального времени; (2) метод вычислительно устойчив в случае плохо обусловленных ковариационных матриц (которые возникают, например, в случае близко расположенных измерений); (3) метод требует задания только относительной ошибки сЦс20 , но не абсолютных значений каждой из ошибок; (4) метод сохраняет основные свойства субоптимальных схем, основанных на фильтре Калмана, таких как оптимальная интерполяция. Последнее означает, что для векторов наблюдений и поля первого приближения /0 и /в величины /, измеренные на одном и том же множестве точек 0<к<К + 1, с К х К ковариационными матрицами: О = ((0 - ? )•(% - ^)), В = ((^в - )•$ - )), 0 <,,у < К , процедура (1)-(2) минимизирует
функционал качества:
1Д = (,+1[| + ОТ ['о-'»], (4)
если параметр сЦсО постоянен во всей области и ковариационная матрица О(.,у) может быть представлена как О(.,у’) = о\™ (.,к), где ^(.,к) определено соотношениями (3). Тот же функционал
(4) минимизируется в фильтре Калмана, который, однако, лишен тех ограничений, которым подвержен представленный метод ИОР - постоянство относительной ошибки и специальный вид ковариационной матрицы.
Полный цикл ассимиляции состоит из двух шагов. На шаге прогноза вычисление поля первого приближения производится путем численного решения конечно-разностных уравнений, аппроксимирующих уравнения переноса концентрации. На шаге ассимиляции значения поля первого приближения обновляются с помощью алгоритма (1)-(3). Следует отметить, что
применимость метода ИОР ограничена сравнительно однородными полями концентрации. В случае, когда поле концентрации характеризуется большими градиентами (например, в течение короткого промежутка времени после начала выброса), применение этого метода может приводить к появлению отрицательных значений концентрации, поэтому на начальных стадиях распространения загрязнения после формирования атмосферных выпадений этот метод неприменим.
3. Результаты расчетов
3.1. Моделирование радиоактивного загрязнения Черного моря в результате Чернобыльской аварии
Для расчета радиоактивного загрязнения Черного моря использовалась карта батиметрии, которая интерполировалась на сетку модели с шагом 27,7 x 23 км и общим числом вычислительных точек 50x30. Вертикальное разрешение модели - 21 уровень в сигма-системе координат. Использовались среднемесячные климатические данные из климатологической базы данных Хеллермана и Розенштейна [10]. Компоненты напряжения поля ветра и температуры интерполировались на расчетную сетку. Использовались климатологические данные о сезонных расходах воды в реках Дунай и Днепр. Пролив Босфор учитывался с помощью гидравлической двухслойной модели, описывающей вток в нижней части моря и выток в верхней части [11]. Поля скорости, температуры, солености и концентрации радионуклидов рассчитывались за период с января 1986 г. по декабрь 1991 г. с временным разрешением 20 минут.
137
На рис. 1 показано распределение Оэ в поверхностном слое в июне 1986 г, рассчитанное на основании данных измерений [10] после выпадений радионуклидов в результате Чернобыльской аварии. Результаты численных экспериментов с использованием ассимиляции данных показаны на рис. 2-3. Радиус влияния полагался равным Rcor = 60км. Наблюдения, обозначенные ромбами,
использовались в процедуре ассимиляции, в то время как наблюдения, обозначенные квадратами, использовались для вычисления среднеквадратической ошибки прогноза.
«О - 600
¿по - ¿во
350 - ¿ПО 300 - 350 250 - 300 200 - 250 150 - 200 100- 150 Э0 - 100 0-50
Рис. 1. Концентрация 137Сэ (Вк/м3) в поверхностном слое после выпадения радионуклидов в результате
Чернобыльской аварии, июнь 1986
27° Е 30° Е 33° Е 36°Е 39°Е 42"Е
к
ГС
ь
Ф Ю I- X
0 О
1
н
о
На рис. 2 показана
зависимость от времени относительной
среднеквадратической ошибки прогноза
*=, К С - с:, )7е( П с:
с
/=1 •• II
использованием и без
использования ассимиляции. Как видно из этого рисунка, ассимиляции на качество полей концентрации
возрастает со временем от начала расчета до ( < 1000 дней. Для больших времен ( < 1000 дней
относительная ошибка поля концентрации радионуклидов б, рассчитанная без ассимиляции,
приблизительно в два раза превышает ошибку поля, рассчитанного с ассимиляцией. Значения
дни
Рис. 2. Относительная среднеквадратическая ошибка прогноза: штриховая линия - без ассимиляции, сплошная влияние линия - с ассимиляцией
рассчитанных
45°Ы-
-45° N
42" М-
-42° N
Л 1 2 з
абсолютной ошибки поля, рассчитанного без ассимиляции б1 =./Х_С -С'оЬ,) , на 10-20 Вк/м3
V /=1 п
превышают значения ошибки поля, рассчитанного с ассимиляцией.
27°Е 30°Е 33°Е 36°Е 39°Е 42°Е 120 - 134
Ш 107 - 120 Ш93 - 107 И ВО-93
■ 67-80 И 53-67
■ 40 - 56!
27-40 13-27 0 - 13
1111 - 124 199 - 111 187-99 174-87 162-74 49 - 62 37-49 25-37 12 - 25 0- 12
45°Ы
42" Ы
Ч5"Ы
-42" N
42" Е
Рис. 3. Концентрация 137Сэ (Бк/м3) в поверхностном слое Черного моря, июнь 1987: а) без ассимиляции; б) с ассимиляцией
у, м
120000-
Рассчитанные без использования и с
10000С
20000
40000
80000
60000
использованием ассимиляции двумерные поля концентрации представлены на рис. 3. На этих рисунках видно резкий максимум концентрации радионуклидов в средней части Черного моря, который уменьшается на другой части рисунка, где представлено распределение концентрации, рассчитанное с ассимиляцией. Местонахождение максимума также смещается под
0
влиянием ассимиляции.
0
40000
80000
120000
Рис. 4. Карта томографии Каховского
водохранилища. Точками обозначено 3.2. Моделирование радиоактивного загрязнения
Метод ассимиляции данных тестировался также на примере сценария потенциального аварийного выброса из Запорожской АЭС и последующего выпадения радионуклидов на поверхность Каховского водохранилища. Расчеты проводились для отрезка времени от 60 по 161 день среднестатистического года. Батиметрическая карта водохранилища интерполировалась на сетку модели с разрешением 1,49x1,1251 км (рис. 4), общее количество узлов модели было равным 101x101. Использовались среднегодовые значения сдвигового напряжения ветра и температуры воздуха, которые интерполировались на сетку модели.
Для тестирования метода ассимиляции данные измерений имитировались результатами расчетов в точках, представленных на рис. 4 для сценария “истинного” выброса. Эти данные затем использовались при проведении расчетов по сценарию “оценочного” выброса для уточнения результатов прогнозирования.
В качестве “истинного” выброса был взят сценарий, который использовался во время плановых учений по ядерному аварийному реагированию, проведенных НАЭК “Энергоатом” в августе 2002 г. Его характеристики были следующие: выброс из Запорожской АЭС 137Сэ в течение 1 часа с интенсивностью 1,4 Е+12 Вк/с. Скорость ветра задавалась постоянной по величине и равной
1.5 м/с, в то время как направление ветра менялось со временем от 120 град. до 150 град. и обратно в течение часа после выброса. Категория устойчивости атмосферы была выбрана нейтральной, осадки отсутствовали. “Истинное” поле атмосферных выпадений на поверхность водохранилища рассчитывалось с помощью модели атмосферного переноса ^Мрирр [2], используемой в системе РОДОС. Временные зависимости концентраций в точках, показанных на рис. 5, записывались с шагом, равным шагу интегрирования модели, и затем использовались как псевдоизмерения в сценарии «оценочного выброса». Этот сценарий характеризовался меньшей в
1.5 раза интенсивностью и тем же временем выброса. Метеорологические параметры были те же, что и для первого сценария. При расчетах с ассимиляцией использовались данные выпадений «оценочного» сценария и данные псевдоизмерений, полученные при расчетах “истинного” выброса.
местонахождение псевдо-измерении,
используемых в ассимиляции
Каховского водохранилища в результате гипотетического выброса на Запорожской АЭС
100000
80000
60000
40000
20000
0
0
40000 х, м
80000
Бк/м**3
I
2800 2600 2500 2450 2400 2000 1 600 1 200 800 400 0
Рис. 5. Концентрация радионуклидов на 75 день после выброса: а) “истинное” распределение; б) без ассимиляции; в) с ассимиляцией(Я_ =6км)
Рис. 6. Концентрация радионуклидов на 125 день после выброса: а) “истинное” распределение; б) без
ассимиляции; в) с ассимиляцией (Ксог = 6км)
Пространственный радиус влияния в процедуре ассимиляции данных полагался равным 6 км, а временной - 3 дня. Отметим, что для случая расчетов распространения радионуклидов в Черном море пространственный радиус влияния был в 10 раз больше. Такое отличие значений этого параметра связано с отличиями горизонтальных геометрических масштабов обоих водоемов. В общем случае значение радиуса влияния должно выбираться (по аналогии с горизонтальным масштабом турбулентных пульсаций) приблизительно пропорциональным характерному геометрическому масштабу водоема.
Результаты расчетов распространения поля загрязнения в водохранилище в случае “истинного” выброса, “оценочного” выброса без ассимиляции и “оценочного” выброса с ассимиляцией представлены на рис. 5-6. Как видно из представленных на рис. 5-6 результатов, процедура ассимиляции приводит к существенному улучшению точности полей концентрации, особенно для более длительных сроков прогноза. Рассчитанные на 125-й день после выброса максимальные концентрации, полученные с использованием ассимиляции, практически совпадают с “истинными”, в то время как максимальные концентрации, рассчитанные без ассимиляции, в 1,5 раза превышают “истинные” значения. Улучшение результатов ассимиляции с увеличением срока прогноза объясняется уменьшением градиентов концентрации.
4. Выводы
В работе представлен метод ассимиляции данных “итерации к оптимальному решению”, используемый в модели распространения радионуклидов в водоемах ТРИТОКС системы РОДОС. Представленный метод, являющийся модификацией известного метода оптимальной интерполяции, прост в реализации и вычислительно эффективен.
Метод верифицирован на основании измерительных данных загрязнения Черного моря в результате радиоактивных выпадений после Чернобыльской аварии. Использование усвоения данных измерений в этом случае привело к существенному уменьшению статистических характеристик ошибок прогноза концентрации 137Сэ в Черном море. Метод ассимиляции данных
137
верифицировался также на примере расчета распространения Сэ после выпадения на поверхность Каховского водохранилища вследствие гипотетического выброса из Запорожской АЭС. В этом случае использование ассимиляции данных также существенно приближало рассчитанное поле концентрации к “истинному”. Метод лучше работает в случае более однородно распределенных полей концентрации.
Из результатов численных экспериментов, представленных в работе, следует, что предложенный метод ассимиляции данных может быть использован для повышения точности прогностических расчетов в моделях водных систем, используемых в системах поддержки реагирования на радиационные аварии.
БЛАГОДАРНОСТИ
Представленная работа выполнялась в рамках проекта 5-й рамочной программы Европейской Комиссии “РДОМЕМ”, контракт № РІКР-СТ-2000-0025. Авторы выражают благодарность В.И. Кошебуцкому за помощь в проведении расчетов, а также В.С. Мадеричу и Д.М. Трибушному за участие в обсуждении результатов.
СПИСОК ЛИТЕРАТУРЫ
1. Margvelashvily N., Maderich V. and Zheleznyak M. THREETOX - A Computer Code to Simulate ThreeDimensional Dispersion of Radionuclides in Stratified Water Bodies // Radiat. Prot. Dosimetry. - 1997. - Vol. 73. -P. 177-180.
2. Ehrhardt J. The RODOS System: Decision Support for Off-Site Emergency Management in Europe // Radiat. Prot. Dosimetry. - 1997. - Vol. 73. - P. 41-44.
3. Zheleznyak M., Heling R., Raskob W. Hydrological dispersion module of the decision support system RODOS // Radioprotection - colloques. - 2002. - Vol. 37. - P. 683-688.
4. Rojas-Palma C., Madsen H., Gering F., Puch R., Turcanu C., Astrup P., Müller H., Richter K., Zheleznyak M., Treebushny D., Kolomeev M., Kamaev D. and Wynn H. Data assimilation in the decision support system RODOS // Radiat. Prot. Dosimetry. - 20O3. - Vol. 104. - P. 31-40.
5. Ahlbrecht M., Borodin R., Borzenko V., Ehrhardt J., French S., Shershakov V., Sohier A., Trakhtengerts E., Verbuggen A. Decision Support Issues in RODOS: The Needs of Decision Makers // Radiat. Prot. Dosimetry. -1997. - Vol. 73. - P. 81-84.
6. Daley R. Atmospheric Data Analysis, Cambridge University Press. - 1999. - 455 p.
7. Robinson A.R., Lermusiaux P.F.J., Sloan N.Q. Data Assimilation / eds. H. Brink, A. Robinson. - John Wiley
&Sons. - 1999. - Vol. 10. - P. 541-593.
8. Kivva S. Data assimilation for 2D model of the radionuclide transport by surface water flow. In: Computational Science - ICCS 2003 / Peter M.A. Sloot, et. al., eds). - Springer. - 2003. - V. II. - P. 619-628.
9. Ковалец И. Ассимиляция метеорологических данных в системах поддержки принятия решений // Математические машины и системы. - 2005. - № 2. - C. 155-167.
10. Veleva B., Koziy L., Yushchenko S., Maderich V., Mungov G. Assessment of radionuclide contamination in the Black Sea using Poseidon / Rodos system // Radioprotection-Colloques. - Vol. 37-C1. - P. 827-832.
11. Мадерич В.С., Никишов В.И., Стеценко А.Г. Динамика внутреннего перемешивания в
стратифицированной среде. - Киев: Наукова думка, 1988. - 240 с.
