CC BY
20
УДК 519.612:632.4
А.В. Гончарский, С.Л. Овчинников2, С.Ю. Романов3 ОБ ОДНОЙ ЗАДАЧЕ ВОЛНОВОЙ ДИАГНОСТИКИ
Работа посвящена разработке методов и алгоритмов решения нелинейной трехмерной обратной задачи волновой диагностики неоднородностей в однородных средах в приближении уравнения Гельмгольца. В силу вычислительной сложности задачи ее решение возможно только с использованием суперЭВМ. Проблема состоит в поиске неизвестного коэффициента в уравнении с частными производными. Задача трансформируется к системе операторных уравнений относительно неизвестного коэффициента и волнового поля. Приведены эффективные алгоритмы ее решения с использованием суперкомпьютеров с параллельной архитектурой.
Ключевые слова: параллельные вычисления, нелинейные волновые модели, обратные задачи.
1. Введение. Томографические исследования (X-RAY, MRT) в настоящее время являются неотъемлемой частью большинства сложных медицинских исследований. Рентгеновская томография широко применяется в технической диагностике. Томографическая аппаратура позволяет получать предельно высокое разрешение.
Однако существует большой класс задач вычислительной диагностики, где томографические методы пока используются чаще на уровне математического моделирования, чем практического применения. К таким задачам относятся прежде всего задачи исследования приповерхностных слоев Земли. Разработанные методы диагностики в настоящее время направлены в основном на исследование слоистых структур, залегающих параллельно или наклонно к поверхности Земли [1-5].
В настоящей работе авторы делают попытку оценить возможности решения обратной задачи диагностики — поиска локальных неоднородностей в однородных, а в перспективе и слоистых средах. Такие методы могут найти применение не только в инженерной геофизике, физике Земли, исследовании структуры вулканов, но и в медицине.
Назовем основные проблемы, возникающие при решении подобных задач:
1) волновая природа источников;
2) необходимость использования суперкомпьютеров;
3) несовершенство используемых моделей.
Остановимся более подробно на каждом из пунктов. Для исследования приповерхностных слоев Земли широко используют волновые источники излучения (акустические волны, сейсмические волны, электромагнитное излучение). Волновые модели, используемые для решения прямых и обратных задач, должны основываться на общепринятых моделях распространения волн в среде. Самая простая модель — скалярное волновое уравнение Гельмгольца [6]. Эта модель учитывает процессы дифракции, рефракции волн, переотражения и т. п. Основной проблемой является то, что, как правило, обратные задачи в таких моделях приводят к нелинейным задачам с очень большим числом неизвестных. Решать их можно только на суперЭВМ [7-10].
В настоящей работе предложены вычислительные алгоритмы решения обратной задачи поиска неоднородности в рамках модели Гельмгольца. Проведенные модельные расчеты показали, что рассматриваемые модели позволяют выявить неоднородность, расположенную в приповерхностных слоях Земли.
К сожалению, реальные задачи намного сложнее. В среде распространяется не одна волна, а несколько (продольные, поперечные, головные, поверхностные и т.д.), т.е. реальная задача является тензорной. В настоящее время в некоторых публикациях сделаны попытки решения прямых задач распространения волн в тензорных средах [11]. Решение обратных задач в тензорной модели не под силу и современным суперЭВМ.
1НИВЦ МГУ, проф., д.ф.-м.н., e-mail: goncharQsrcc.msu.ru
2НИВЦ МГУ, электроник, e-mail: oslQstarnet.ru
3НИВЦ МГУ, с.н.с., к.ф.-м.н., e-mail: romanovQsrcc.msu.ru
В некоторых случаях благодаря тому, что скорости разных волн в среде различны, их можно разделить по времени прихода. Если это удается сделать, то для каждой отдельной волны можно использовать скалярную модель. Именно в этой модели — скалярной волновой модели Гельмгольца — мы и будем продолжать рассмотрение.
2. Математическая постановка обратной задачи. Рассмотрим задачу, описываемую следующим скалярным волновым уравнением:
Au(r,q,w) + n2(r,w)u(r,q,w) = f(r,q,w), (1)
где гей3 — точка в пространстве, А — оператор Лапласа по переменной г. Уравнение (1) в скалярном приближении описывает акустическое или электромагнитное поле u(r, д, ш), генерируемое источником, описываемым функцией /(г, д, ш). Будем предполагать, что источник является точечным гармоническим осциллятором, расположенным в точке q € Q. Тогда эта функция имеет вид /(г, д, ш) = — ш — круговая частота источника. Пусть неоднородность среды вызвана только изменениями фазовой скорости с(г), имеет место соотношение к(г,ш) = ш/с{г). Вне области неоднородности будем предполагать для простоты к(г,ш) = kq = ш/cq, где со = const известна. Рассмотрим функцию Грина
G(r,q,u) = -—^-г-ехр (t— \\r - д|| ) , (2)
4-7Г \\г — g|| V Со J
2
которая удовлетворяет уравнению AG(r, g, ш) + g, ш) = —S(r — q).
co
Предположим, что поле u(r,q,u) измеряется в области /'. г (г /'. т.е. приемники излучения "пробегают" область Р. Точечные источники излучения "пробегают" область Q, q G Q, неоднородность расположена в ограниченной области R. Раздельно записывая уравнения для областей ДиР,с учетом (2) получаем нелинейную систему уравнений [12, 13]
{u(r, д, ш) = щ(г, д, ш) + ш2 f G(r', г, ш)£(г')и(г', д, ш) dr', г € R,
R (з)
U(p,q,u) = ш2 f G{r',р,ш)£(г')и(г',д, ш) dr', р € Р,
я
где введены обозначения U(p,q,u) = и(р, д, ш) — щ(р, д, ш),
uo(r,q,u) = - J G(r', г, ц?)/(г', д, ш) dr', ^(r) = с^2 - с~2(г).
Q
В обратной задаче неизвестными считаются свойства среды £(г) и поле u(r,q,u). Известными считаются уравнения связи (3), в которых функция U(р, д, ш) получается как результат измерений в области расположения приемников. Требуется по этим измерениям восстановить неоднородность среды £(г).
3. Алгоритмы решения обратной задачи. В основе итерационного метода решения обратной задачи лежит регуляризованная процедура Гаусса-Ньютона [12, 13]:
арЕ)~1(Р'*Р(Хр)
Хр+1 — Хр
t{f';f'p
Ор(Хр - С)),
(4)
где Хр = (щ ,... , Ср) — вектор неизвестных, I = q ■ d, q — число используемых частот зондирования, й — число положений источников, Е — единичный оператор, 0 < г < 1, ар — параметр регуляризации, С — вектор регуляризации, ¿^ = Р'(Хр) — производная Фреше, Рр — оператор, сопряженный р — номер итерации. Формулы для Р(Хр), ^ = Р'(Хр) и РрРр приведены в [8]. В матрично-векторном виде матрица Рр имеет размер (п + т)1 х п(1 + 1), эрмитова положительно-определенная матрица РрР'р размера п{Ь+ 1) хп(|+ 1) имеет вид
F'*P' =
р Р
'С1*С1 О
о
DUC1
О
с2*с2
о
D2*C2
о о
си с*
БиСг
Di*Di
Cl*D1 C2*D2
ct*Dt D2*D2 +
D**D*
где п — количество точек сетки в трехмерной расчетной области Д, содержащей неоднородность; т — набор точек измерения добавочного поля II; С и Б — некоторые матрицы [8]. Начальное приближение
V ( (!) 0) п\ 0)
возьмем в виде Ха = (щ ,... , и), где иа — поле, полученное при заданных значении частоты и положении источника.
Наиболее вычислительно трудоемкой операцией в рассматриваемой процедуре Гаусса-Ньютона (4) является вычисление матрицы и обращение уравнения Сх = Ь на каждом шаге итераций, где С = (РрРр + арЕ), Ь = Рр*Р(Хр) + ар(Хр — С), С = С* > 0. Для сравнения алгоритмов нами были разработаны однопроцессорные варианты программы обращения оператора {Р'^Р' + арЕ) методами Гаусса, Холецкого, простой итерации, нестационарным трехслойным методом и методом сопряженных градиентов [14, 15]. Далее для алгоритма, показавшего наибольшую эффективность, был разработан многопроцессорный вариант программы с использованием технологии МР1.
Асимптотика алгоритмической сложности метода сопряженных градиентов имеет вид
• 0((п3)2), где — число итераций, п — число точек области решения по одной координате. Для методов простой итерации и нестационарного трехслойного метода формула алгоритмической сложности выглядит так же, но Нц — другое число итераций. Таким образом, скорость итерационных процессов зависит от числа итераций. Обычно итерации заканчивают, когда невязка становится равной ошибке входных данных 8. Методы Гаусса и Холецкого являются прямыми методами решения системы линейных уравнений. Их алгоритмическая сложность имеет вид 0((п3)3).
В настоящей работе были проведены численные эксперименты сравнения времени работы пяти рассмотренных алгоритмов на однопроцессорных компьютерах. Рассматривалась задача с двумя источниками на двух частотах, исходные данные были известны на сетке 13x13 точек. Проводилось 5 итераций метода Гаусса-Ньютона (т.е. 5 раз решалась линейная система). Задача решалась на двух сетках 7x7x7 точек (размер матрицы 1715x1715) и 9x9x9 точек (размер матрицы 3645x3645). Итерации заканчивали при достижении невязки 8 = 3• Ю-4 и 8 = 1 • Ю-4. Время расчета (в секундах) для каждого варианта приведено в таблице. К этому времени необходимо добавить время пятикратного вычисления матрицы, которое составило 255 с для матрицы 1715x1715 и 1625 с для матрицы 3645x3645.
Сравнение времени работы различных алгоритмов
Метод Размер матрицы
1715x1715 3645x3645
<5 = 3- Ю-4 <5 = 1- Ю-4 ¿ = 3- Ю-4 <5 = 1- Ю-4
Простая итерация 1113 3566 16102 27697
Нестационарный трехслойный 52 128 224 603
Сопряженных градиентов 29 64 137 304
Гаусса 115 1088
Холецкого 41 295
Как видно из таблицы, метод сопряженных градиентов сходится значительно быстрее метода простой итерации и примерно в два раза быстрее нестационарного трехслойного метода. Прямой метод Холецкого быстрее метода Гаусса. Однако метод Холецкого с увеличением размерности начинает уступать методу сопряженных градиентов. Кроме того, метод сопряженных градиентов, в отличие от прямых методов, обладает очень высоким потенциалом распараллеливания вычислений. Поэтому для дальнейших расчетов на параллельном компьютере был выбран метод сопряженных градиентов.
4. Использование специального вида матрицы. Как было отмечено выше, количество эле-
ментов матрицы решаемой системы уравнений растет примерно как 0(п8). Используя разреженность
и эрмитовость матрицы ¿^ (см. (5)), удалось сформировать уплотненную матрицу с асимптотикой
количества элементов 0(п7). Пусть N = • Ыу • тогда очевидно, что для всех строк матри-
цы меньше йцЫ ненулевыми являются только 2/У элементов. Кроме того, матрица К^К^ является
эрмитовой. Все это позволяет хранить только 2(<1д + 1) "клеток" размером ЖхЖ:
р р
С^С1
с2*с2
сис1
писг
С1* Б1 С2*В2
сивг В2*В2 +
вивг
(6)
Таким образом, при Ж = 213,с^ = 5ид = 30 необходимый объем памяти для хранения матрицы составляет 193 Гб вместо 14 570 Гб для неуплотненной матрицы.
Хранение уплотненной матрицы ведет к усложнению алгоритма решения системы линейных уравнений, однако практически не сказывается на времени решения системы.
5. Организация параллельных вычислений. Как было показано выше, численное решение обратной задачи (3) сводится к итерационному процессу (4), на каждом этапе которого необходимо решить систему линейных уравнений с матрицей (¿^¿^ + арЕ). Основными вычислительными трудностями данной задачи являются большой объем требуемой оперативной памяти для хранения матрицы (¿^¿^ + арЕ) и большое количество собственно вычислительных операций. Как видно из (5), процедура вычисления матрицы Рр*Рр сводится в конечном итоге к перемножению матриц размера ЖхЖ, N = Ых ■ Ыу ■ в количестве Мд раз и составляет порядка 4 • й ■ д • 0(Ж3) операций. Сложность решения системы линейных уравнений также составляет порядка ЩгО(Ы2), где — число итераций.
Одним из возможных путей решения вышеуказанных трудностей может являться использование вычислительных комплексов типа "кластер". При данном подходе возможно как распределить необходимую оперативную память по различным вычислительным узлам, так и организовать на этих узлах параллельные вычисления.
Для решения обратной задачи нами была выбрана система организации параллельных вычислений типа "один ведущий, остальные ведомые". Такая система позволяет сосредоточить всю логику работы программы, чтение/запись данных, выдачу диагностики и сообщений в одном месте — на одном вычислительном узле.
Выделение одного вычислительного узла в качестве "ведущего" позволяет сильно упростить отладку и сопровождение программы. Все логические решения принимаются только на "ведущем" узле, остальные узлы работают в режиме "управления потоком данных". Они ожидают прихода исходных данных для расчета матрицы, после приема этих данных производят расчет своей части матрицы и правой части уравнения и далее вызывают процедуру решения системы линейных уравнений.
"Ведущий" узел вычисляет начальные данные для итераций, читает промежуточные данные из файла сохранения, создает этот файл в случае необходимости, осуществляет вычисления после решения системы линейных уравнений и принимает все логические решения: об окончании итераций, об окончании выполнения программы и т. п.
Ниже описываются методы оптимизации использования оперативной памяти и сокращения времени расчетов.
5.1. Параллельное вычисление матрицы. Основным "потребителем" оперативной памяти является массив матрицы Остальные массивы значительно меньше и могут иметь копии на
каждом вычислительном узле. Например, при тех же параметрах N = 213, й = 5 и д = 30, размер массива правой части системы линейных уравнений и массива решения составляет около 11 Мб (для сравнения: размер матрицы 193 Гб).
Структура используемого алгоритма позволяет разбить массив матриц на части и распределить их вычисление и хранение по отдельным узлам вычислительной системы. Естественным методом распределения массива матрицы по узлам является выделение для каждого узла своего диапазона строк матрицы. В данной реализации принято решение о выделении каждому узлу равного числа строк матрицы. Этот метод не является оптимальным с точки зрения равномерного распределения нагрузки на вычислительные узлы, но сильно упрощает алгоритм синхронизации работы узлов.
В результате данного подхода минимальное количество узлов вычислительной системы, требуемых для решения конкретной задачи, приблизительно равно размеру матрицы в байтах, деленному
на 1 Гб (при оперативной памяти узла, равной 1 Гб). Данное число может быть немного больше или меньше, в зависимости от распределения задачи по конкретным вычислительным узлам.
Матрица формируется сразу в уплотненном виде. Поскольку она является эрмитовой, то вычисление нижней "полосы" матрицы (последние N строк) не производится.
В процедуре расчета матрицы также предусмотрена возможность добавления регуляризирующей диагональной матрицы аЕ.
5.2. Параллельное решение системы линейных уравнений. Для решения системы линейных уравнений был выбран метод сопряженных градиентов. Как видим, основной операцией в данном методе является умножение строки матрицы на вектор, поэтому естественным методом распределения вычислений по узлам кластера является выделение для каждого узла своего диапазона строк матрицы. При такой организации вычислений основной объем итерационных вычислений выполняется независимо на отдельных узлах. Для неуплотненной матрицы реализация данного алгоритма весьма проста, для уплотненной матрицы (6) более сложная. Рассмотрим параллельную работу алгоритма более подробно.
Алгоритм работает с уплотненной матрицей. Сначала проводится чтение промежуточных значений из файла сохранения. Эту операцию выполняет узел с номером 0 и затем широковещательным обменом рассылает прочитанные значения на все узлы. Рассылка осуществляется процедурой МР1_Всаз1;, также данная процедура обеспечивает синхронизацию, т. е. узел ожидает получения данных и не продолжает работу, пока обмен данными не будет завершен.
Далее реализуется собственно итерационный процесс решения системы линейных уравнений. Поскольку каждый узел имеет значения всех элементов массивов Хк и Хк~1, то в каждом узле могут быть независимо вычислены значения СХк для тех строк матрицы, с которыми этот узел работает, за исключением последних N строк. При этом в силу структуры уплотненной матрицы вычисления для Сп*С'\ Сп*П'1 и Вп*Па приходится выполнять отдельно. Поскольку уплотненная матрица не содержит всех элементов для последних N строк, значения СХк вычисляются отдельно с использованием того факта, что матрица Ер*Ер является эрмитовой. Частичные суммы для последних N строк матрицы вычисляются в отдельном массиве в каждом узле. Далее эти частичные суммы передаются в последний узел, где суммируются (с помощью процедуры МР1_Ыес1исе). Вычисленные окончательные суммы рассылаются всем узлам с помощью процедуры МР1_Всаз1;. После вычисления всех элементов СХк каждый узел имеет возможность вычислить Хк+1 для тех строк матрицы, которые содержатся в данном узле. После чего узлы обмениваются вычисленными значениями Хк+1 с помощью процедуры МРЦЗсаБ!;. Очередная итерация завершена.
Важным элементом итерационного алгоритма является критерий прекращения итераций. Нами был выбран критерий, основанный на величине невязки ЦСХ^"1"1 — Ь||. Вычисление невязки производится аналогично вычислению СХк.
Несмотря на высокое быстродействие кластерного суперкомпьютера, из-за большого объема задачи необходимо иметь возможность осуществлять прерывания и продолжать расчет задачи с прерванного места. В связи с этим была реализована процедура сохранения промежуточных результатов. Мы можем сохранять значения вектора X = (г^1-5,..., £), вычисленные на каждом шаге итерации Гаусса-Ньютона. Для того чтобы не испортить вычисленные на предыдущем шаге значения при прерывании задачи в процессе записи, массив сохранения создается двойного размера: сначала сохраняемые значения записываем в первую половину, затем во вторую, затем снова в первую и т. д.
При каждом старте программа проверяет, существует ли файл сохранения: если существует, то промежуточные данные загружаются из файла и счет продолжается с прерванного места, если файл не существует, то он создается и в него записываются исходные значения параметров счета.
Для внутренних итераций решения системы линейных уравнений также используются файлы сохранения. Относительно малый объем промежуточных данных позволяет сохранять эти данные на каждом шаге итерации. Алгоритм сохранения/восстановления промежуточных данных аналогичен алгоритму для внешних итераций.
6. Математическое моделирование. Численные эксперименты были проведены на суперкомпьютере СКИФ МГУ-"Чебышев". Рассматривалась следующая модель. На глубине 45 м от поверхности Земли располагалось два одинаково отражающих объекта на расстоянии 90 м между собой и один такой же объект на глубине 135 м (рис. 1). Фазовая скорость в среде равна со = 1800 м/с. На
поверхности Земли над объектами симметрично в квадрате размером 380x380 м располагалось 40x40 приемников в узлах решетки с равномерным шагом. Источник располагался на поверхности посередине между объектами. Зондирование проводилось на частотах 40 и 60 Гц, что соответствует длинам волн 45 и 30 м. При решении обратной задачи область расчетов R имела размеры 180x180x180 м, в ней была введена сетка 33x33x33 точек. Расчет линейной системы заканчивался при достижении невязки 5 = 3- Ю-4. Время расчета одной внешней итерации на 512 процессорах составило около 7.5 ч, было реализовано 80 итераций метода сопряженных градиентов.
В этом численном эксперименте предполагалось промоделировать возможность решения задачи на большой сетке в области решения. Число неизвестных составило около 110 000 комплексных чисел.
Поверхность
Рис. 1. Исходное изображение трех объектов: сечение в вертикальной плоскости (а); сечение в горизонтальной плоскости, проходящее через два объекта (б); сечение в горизонтальной плоскости, проходящее
через один объект (в)
Поверхность
Рис. 2. Восстановленное изображение трех объектов: сечение в вертикальной плоскости (а); сечение в горизонтальной плоскости, проходящее через два объекта (б); сечение в горизонтальной плоскости,
проходящее через один объект (в)
На рис. 1, а, б, в приведено изображение амплитуды исходных данных в вертикальном и горизонтальных сечениях, проходящих через объекты. На рис. 2, а, б, в — восстановленное изображение в тех же сечениях. Проведенные модельные расчеты показывают, что в рамках поставленной обратной задачи волновой томографии удается восстановить локализованные неоднородности как в непосредственной близости к поверхности, так и на глубине.
7. Заключение. По результатам работы можно сделать следующие выводы.
1. Рассматриваемая обратная задача с вычислительной точки зрения является очень трудоемкой. Количество неизвестных задачи растет примерно как 0(п4), где п — число точек в одном направлении трехмерной сетки. Единственный путь решения этой проблемы состоит в использовании суперкомпьютеров. Использование суперкомпьютеров семейства СКИФ позволяет довести количество вычислительных узлов до 500, в несколько раз поднять размерность рассматриваемой задачи (п > 30).
С увеличением размерности удается рассматривать все более содержательные задачи. Таким образом, использование современных компьютеров позволяет решать обратные задачи волновой томографии в нелинейной постановке коэффициентных задач уравнения Гельмгольца.
2. С точки зрения практических применений основные проблемы связаны с несовершенством используемой модели. Несмотря на то что приближение Гельмгольца описывает эффекты рефракции, дифракции и т.д., скалярная волновая модель далека от действительности, например, в инженерной сейсмике. Все дело в том, что реально мы вынуждены иметь дело не с одной волной, а с целым набором волн (поперечная, продольная, поверхностная, головная и т.д.), которые реально распространяются в среде.
3. Использование тензорных моделей перспективно для решения прямых задач, но для обратных задач нереально даже для современных суперЭВМ.
4. Авторы статьи видят единственную возможность в получении прорывных результатов в поиске таких задач, где физически удается разделить продольные, поперечные, поверхностные волны, используя то обстоятельство, что они имеют разные скорости распространения. При достаточно большой глубине залегания неоднородностей это удается сделать, например, в задаче диагностики внутренней структуры вулканов [16].
СПИСОК ЛИТЕРАТУРЫ
1. Claerbout J. Imaging the Earth's interior. Boston: Blackwell Scientific Publishers, 1986.
2. Robinson E., Durrani T.S., Peardon L. G. Geophysical signal processing. Englewood Cliffs: Prentice Hall, 1986.
3. Gazdag J., Squazzero P. Migration of seismic data // IEEE. 1984. N 72. P. 1302-1315.
4. Yilmaz O. Seismic Data Processing. Tulsa: Society of Exploration Geophysicists, 1987.
5. В ay sal E., Kosloff D. D., Sherwood J. W. C. Reverse time migration // Geophysics. 1983. N 48. P. 15141524.
6. Ramm A. G. Multidimensional Inverse Scattering Problems. London: Longman Group, 1992.
7. Гончарский А.В., Романов С. Ю., X ар чен к о С. А. Обратные задачи диагностики с синтезированной апертурой // Труды Международной научной конференции "Параллельные вычислительные технологии" (ПаВТ'2007). Т. 1. Челябинск: Изд-во ЮУрГУ, 2007. С. 109-114.
8. Гончарский A.B., Романов С. Ю., Харченко С. А. Обратная задача акустической диагностики трехмерных сред // Вычислительные методы и программирование. 2006. 7. № 1. С. 117-125.
9. Виноградов Н. С., Гончарский А.В., Овчинников С. JL, Романов С.Ю. Об одной задаче волновой томографии // Труды Всероссийской научной конференции "Научный сервис в сети ИНТЕРНЕТ: многоядерный компьютерный мир — 15 лет РФФИ". М.: Изд-во МГУ, 2007. С. 129-133.
10. Агаян Г.М., Виноградов Н. С., Гончарский А. В. и др. Диагностика трехмерных сред методами синтезированной апертуры // Вычислительные методы и программирование. 2007. 8. № 2. С. 5-10.
11. Кучунова Е.В., Садовский В.М. Численное исследование распространения сейсмических волн в блочных средах на многопроцессорных вычислительных системах // Труды Международной научной конференции "Параллельные вычислительные технологии" (ПаВТ'2008). Челябинск: Изд-во ЮУрГУ, 2008. С. 130-141.
12. Bakushinsky A.B., Goncharsky A.V. Ill-posed Problems. Theory and Applications. Dordrect: Kluwer Academic Publ., 1994.
13. Bakushinsky А. В., Goncharsky A. V., Romanov S.Yu., Seatzu S. On the identification of velocity in seismic and in acoustic sounding. Firenze, 1994.
14. Рябенький B.C. Введение в вычислительную математику. М.: Физматлит, 1994.
15. Тыртышников Е. Е. Методы численного анализа. М.: Академия, 2007.
16. Глинский Б.М., Собисевич А. Л., Хайретдинов М. С. Опыт вибросейсмического зондирования сложно построенных геологических структур (на примере грязевого вулкана Шуго) // Докл. РАН. 2007. 413. № 3. С. 1-5.
Поступила в редакцию 29.04.09
