Численная реконструкция трехмерной скорости звука методом граничного управления Текст научной статьи по специальности «Математика»

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

______________________________2011 г. Выпуск 3 (22). С. 92-98___________________________

УДК 517.9

ЧИСЛЕННАЯ РЕКОНСТРУКЦИЯ ТРЕХМЕРНОЙ СКОРОСТИ ЗВУКА МЕТОДОМ ГРАНИЧНОГО УПРАВЛЕНИЯ

Л. Н. Пестов, В. М. Болгова, А. Н. Данилин

Введение

Пусть О - ограниченная область в К” (п > 2) с гладкой границей Г. Задача, которую мы будем называть прямой, - начально-краевая задача для волнового уравнения:

А-ыи -Аы = 0 тОх(0,Т), (1)

с

ы\г=0= Ы '=0=0 шО, (2)

ы \гх0,Т ]" = I, (3)

где с(х) - положительная функция (скорость звука); (-)у - нормальная производная; / -граничное управление (Неймана); и = и1 (х, ') - решение (волна). В силу конечной скорости распространения волн в любой момент времени имеет место включение

Бирры1 (•, ') сО', ' >0,

где

О5 ={х еО \ё1Б1;Дх,Г) < 5'}, я > 0,

и расстояние понимается как время распространения волны. Подобласть О5 - суть часть О, заполненная волнами к моменту времени ' = 5. В частности, для Т > Т* = Бир{ё181;Дх, Г) \ х е О},

имеет место равенство ОТ = О. С системой (1)-(3) связан оператор реакции ЯТ, который действует по правилу

ЯТ I = и1 \

^ J и 1Гх[0,Т] •

Рассмотрим систему (1)-(3) с финальным моментом ' = 2Т ; пусть Я2Т - соответствующий оператор реакции. В силу принципа причинности, следующего из конечности скорости распространения волн, оператор Я2Т определяется значениями скорости звука в подобласти ОТ и не содержит информации о с |о ^ . Это мотивирует следующую постановку обратной

динамической задачи: для фиксированного момента Т > Т* по заданному оператору Я2Т найти скорость звука во всей области О .

Мы используем метод граничного управления (ВС-метод), предложенный М. И. Бели-шевым [1, 2] и его версию, предложенную в [3]. Настоящая работа является продолжением работы [4], где была развита версия [3] и численно протестирована в двумерном случае. Мы представляем результаты численного моделирования трехмерной обратной динамической задачи. Мы также формулируем основные математические утверждения нашего подхода и отсылаем к статье [4] за деталями и доказательствами. Всюду в статье предполагается выполненным неравенство Т > Т*.

Энергетические формы и задача граничного управления

Билинейные формы, которые мы вводим в этом разделе - основной инструмент решения обратной задачи. Волна Ы зависит от управления / линейно. Тогда кинетическая и потен-

циальная энергия в момент времени ' = Т - квадратичные функционалы в пространстве управлений. Введем соответствующие билинейные формы:

к (/, &) = |О 4) и/(x, Т (x, Т )№х;

'О с2( х)

Р(/, &) = \О(У и1 (х, Т), Уы&(х, Т))Лх. (4)

Вместо кинетической формы будем использовать следующую энергетическую форму:

С (/, &) = {О -зт-ты1 (х,Т )и& (х,Т )№х, (5)

•,О с (х)

которую называем скалярным произведением финальных состояний (следуя теории управлений, мы называем волну в финальный момент и1 (., Т) финальным состоянием). Замечательный факт состоит в том, что все эти формы явно определяются данными обратной задачи (т. е. оператором реакции Я2Т ):

С (/ ’ &)=Цт-/(Я 2Т&)+!/ - 2 & (ш2Т/)]й^Г> Р(/, & )=.и ][. г | (Я 2Т&)++2 & | (Я 2Т/)]<МГ,

где

а±(,') = [«(•,') + а(,2Т - ')]/2,

(!а)(•,') = .а(•,я)<&, ' е 0,2Т].

Другой ингредиент ВС-метода - приближенная граничная управляемость динамической системы (1)—(3). Это означает, что произвольная функция может быть сколь угодно точно приближена некоторым финальным состоянием и/ (., Т) в некоторой естественной норме (т. е. использующей конечное число производных). Другими словами, уравнение (задача граничного управления)

и/(х,Т) = р( х) (6)

относительно управления / плотно разрешимо для любой функции р в некоторых естественных функциональных пространствах (скажем, £2(О)). Вообще говоря, для решения задачи граничного управления (6) необходимо знать скорость звука. Но другой замечательный факт (помимо отмеченных выше явных представлений (4), (5)) состоит в том, что для любой гармонической р (т. е. Ар = 0 ) приближенное решение уравнения (6) может быть найдено только при заданном операторе реакции и неизвестной скорости звука. Более того, как показано в [4] можно контролировать близость финального состояния к заданной гармонической функции «не сходя» с границы. В [4] был предложен следующий алгоритм реконструкции скорости:

1. Задача граничного управления для гармонических функций. Фиксируем гармоническую функцию р и находим соответствующее управление /р, такое, что имеет место (приближенно) (6). Найти такое управление можно путем минимизации квадратичного функционала

Ф(/) + ' "2

где

(Я2Т/)(-,Т)-р

£2(Г)

Фактически это линейная задача. Заметим, что функционал Ф определяется данными обратной задачи в силу (4):

Ф(./) = !О\^и/(.,Т)-Ур\2 йх.

Решаем эту задачу для некоторого множества гармонических функций.

2. Реконструкция р = 1/ с2. Пусть р и у - произвольная пара гармонических функций и /р и /у - соответствующая пара управлений из предыдущего пункта. Тогда мы имеем равенство (см. (5))

|ОР(х)р(х/(х)йх = С(/р, &/) . (7)

Известно, что множество всех произведений гармонических функций плотно в Ь2 (О). Следовательно, функция р единственным образом определяется из (7).

Заметим, что рассматриваемая обратная задача нелинейна. Тем не менее, предложенный алгоритм использует только две линейные процедуры. Как следствие, не возникает проблем с минимизацией невыпуклых функционалов, что часто предлагается при оптимизационных подходах к решению нелинейных обратных задач.

Ниже рассматривается аппроксимация Галеркина для решения прямой задачи и «проекция» обратной задачи на конечномерное пространство.

Аппроксимация Галеркина и дискретная обратная задача

Пусть О - ограниченная область в К3 с триангуляционной сеткой, содержащей N вершин и NА тетраэдров. Пусть пространство граничных управлений Nc -мерно с базисными

управлениями /_,..../N . N кусочно-линейных базисных функций у,...,/ на О определены условиями у(х}.) = 5^, где 5^ - символ Кронекера, х1,...,xN - узлы (вершины) сетки. Рассмотрим аппроксимацию Галеркина решения прямой задачи:

и.N(х')=Е != иа (')/”(х), и(0) = 0, и'(0) = 0,

Тогда матрица и (') - суть решение задачи Коши для системы обыкновенных дифференциальных уравнений с постоянными коэффициентами:

ми"(') + ки = ^ ('), и(0) = и'(0) = 0,

где М и К - матрица масс и матрица жесткости:

м = р(х)у(х)у,-(x)dx, К = (V у(хХ V у^(х))йх,

а ^(') - матрица-функция управлений (размера N х Nc),

Ъа(* ) = |Г у, ( х) /а( х ' )й Г .

Заметим, что для каждого узла х}. имеет место равенство:

и/( х,')=иа(').

Мы используем кусочно-постоянную модель скорости: с(х) = ск в к -ом тетраэдре. Сформулируем дискретную обратную задачу: по заданной матрице-функции

иа('), I е Е,' е [0,2Т], (8)

(где Е - множество граничных узлов) найти значения ск,к = 1,...,NА .

Пестов Л. Н., Болгова В. М., Данилин А. Н. Численная реконструкция трехмерной скорости звука... Главный инструмент решения этой задачи в нашем подходе - следующие матрицы:

C = U*(T)MU(T), P = U*(T)KU(T),

где звездочка означает транспонирование. Эти матрицы - дискретные аналоги энергетических билинейных форм (3), (4) и также явно выражаются по данным (8) обратной задачи. Ранг матрицы C «измеряет» (граничную) управляемость, поскольку rankC = rankU(T). Равенство rank U = N означает управляемость и имеет место только при достаточно большом количестве управлений. Таким образом, по рангу матрицы C можно судить о том, управляема ли система на данной сетке и при имеющемся количестве управлений или нет (заметим, что в дискретной постановке нет смысла говорить о приближенной управляемости). Алгоритм решения обратной задачи в дискретной постановке состоит в следующем:

1. Задача управления. Естественным аналогом гармонической функции на сетке является вектор Y, удовлетворяющий равенству K Y = 0. Используя матрицы C и P, можно найти матрицу-функцию Fharrm (t) управлений [4], которая порождает известные линейно-независимые гармонические состояния Фа, а = 1,..., Nb, где Nb - количество граничных узлов. Решение задачи Коши:

MV(t) + KV = Fharm (t), V(0) = V'0) = 0 обладает ключевым свойством:

V (T) = Ф,

где Ф - известная матрица размера NхNb, столбцы которой Фа,а = 1,...,Nb линейно независимы и принадлежат ядру матрицы жесткости K, KФа = 0 .

2. Реконструкция скорости. В результате решения задачи управления приходим к равенству [4]:

Ф*М Ф = Ch , (9)

harm ’ v '

где матрица Char^m определяется по Fharrm и C. Уравнение (9) - суть система линейных алгебраических уравнений относительно значений рк = 1/ck2, к = 1,...,NA , поскольку

где Ак - к -ый тетраэдр. Заметим, что матрица системы не зависит от данных обратной задачи, а полностью определяется только выбранной сеткой. Данные обратной задачи определяют только правую часть (9).

Численное моделирование

В качестве области О был выбран цилиндр с радиусом 10 см и высотой 20 см. Количество тетраэдров - 37620. Время наблюдения - 2Т, Т = 0,15 мс. Граничные трансдьюсеры (являются одновременно источниками и приемниками) в количестве 210 равномерно расположены только на боковой поверхности цилиндра (рис. 1).

Рис. 1. Сетка и расположение трансдьюсеров

Модель скорости представлена на рис. 2 (1000т / с < с( х) < 3000т / с ). Заметим, что в модели имеет место общее изменение скорости вдоль вертикального направления.

Рис. 2. Модель скорости

Базисные управления выбирались как граничные точечные источники с задержанным по времени импульсом Рикера:

/ав(Х О = 5(x, Xa)r (t-МХ в = I...100,

где 5(х, ха) - граничная дельта-функция, сосредоточенная в граничном узле ха, А - шаг временной задержки, r(t) - импульс Рикера с доминирующей частотой f = 60MHz (рис. 3).

Рис. 3. Импульс Рикера

Результат реконструкции скорости представлен на рис. 4. Относительная погрешность реконструкции в евклидовой норме составляет 7 %. Использовалось различное число гармонических функций. Чем больше это число, тем лучше реконструкция. Рис. 4 соответствует 370 гармоническим функциям, что дает 370 х 371/2 = 68635 уравнений относительно 37620 неизвестных в (9). Тем не менее, матрица системы оказывается вырожденной и при решении системы (9) проводилась регуляризация.

Рис. 4. Реконструкция скорости Заключение

В работе описана первая попытка решения трехмерной обратной динамической задачи с помощью метода граничного управления. Основные выводы следующие: 1) линейных характер BC-решения дает возможность эффективно использовать библиотеки линейной алгебры, типа LAPACK; 2) если время наблюдения достаточно велико, задача граничного управления для гармонических функций численно корректна и удовлетворительная точность ее решения достигается при разумном числе граничных управлений; 3) задача реконструкции скорости требует большого числа гармонических функций (для того, чтобы ранг системы (9) равнялся числу неизвестных Nд ). Непроверенной остается гипотеза о том, что ранг Nд достигается,

если число гармонических функций взять равным числу граничных узлов Nъ (это максимально возможное количество гармонических функций).

Увеличение числа гармонических функций приводит к очень большим системам уравнений с плотными матрицами и требования на вычислительные ресурсы резко возрастают. Выход видится в использовании более грубых конечно-постоянных моделей скорости (увеличение размера ячеек, в которых скорость постоянна), что уменьшает размер системы.

ЛИТЕРАТУРА

1. Belishev, M. I. Recent progress in the boundary control method. Inverse Problems, - Vol. 23. -5. - pp. 1-67. - 2007.

2. Belishev, M. I. Boundary Control Method in Dynamical Inverse Problems. An Introductory Course by M.I.Belishev. Gladwell G.M.L., Morassi A., Editors . Dynamical Inverse Problems: Theory and Application. CISM Courses and Lectures. - Vol. 529. - Wien, Springer. - p. 85150. - 2011.

3. Pestov, L. N. On reconstruction of the speed of sound from a part of boundary. Journal of inverse and ill-posed problems. - 1999. - Vol. 7. - 5. - pp. 481-486. - 1999.

4. Pestov, L., Bolgova, V. and Kazarina, O. Numerical recovering of a density by the BC-method. Inverse Problems and Imaging. - Vol. 4. - 4. - pp. 701-712. - 2010.