Журнал вычислительной математики и математической физики, 2020, T. 60, № 4, стр. 652-662

1. ВВЕДЕНИЕ

Развитие методов численного моделирования акустических волновых полей мотивируется их приложениями в задачах обращения (см. [1]–[4]). Обратные задачи часто формулируются как поиск минимума функционала, который замеряет уклонение значений смоделированного волнового поля от зарегистрированных данных в приемниках. “Смоделированное волновое поле” означает результат расчета прямой задачи для заданной модели (скорости звука в среде). Модель является основным аргументом функционала, и способ решения обратной задачи по сути представляет собой отыскание значения аргумента, доставляющего функционалу минимальное значение. Эффективность алгоритма обращения напрямую зависит от производительности и эффективности алгоритмов решения прямой задачи, которое требуется на каждом шаге итерационного процесса минимизации.

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

Другая сторона того же самого вопроса – это точность получаемого решения прямой задачи. В качестве соответствующей метрики мы используем величину дисперсионной ошибки. Эта метрика, вычисляемая для случая однородной среды, с определенными предостережениями применима и для неоднородных сред, если параметр ppw (число точек на длину волны) рассматривается в определенном интервале значений. Тривиальным способом уменьшения дисперсионной ошибки является использование более мелких шагов разностной сетки, но такой путь может оказаться неприемлемым в силу непомерных затрат.

Для уменьшения дисперсионной ошибки мы используем подход, описанный в [5] (см. также [6]). Трехмерное уравнение Гельмгольца аппроксимируется на 27-точечном шаблоне, полученном путем комбинации нескольких стандартных шаблонов. Существенно, что аппроксимация значения функции вне оператора дифференцирования также происходит на всем шаблоне. Таким образом, получается пяти-параметрическое семейство разностных схем. Значения параметров можно подобрать так, что в интересном для приложений диапазоне значений ppw дисперсионная ошибка полученной схемы второго порядка оказывается значительно лучше, чем для стандартной схемы. Полученная таким образом оптимизированная разностная схема по своей структуре и свойствам близка к схемам повышенного порядка аппроксимации (см. [7]–[10]), но наш подход выгодно отличается своей регулярностью – коэффициенты получаются в результате работы алгоритма оптимизации.

Для решения задачи оптимизации мы комбинируем классический градиентный метод и симплекс-метод (см. [11], [12]), не требующий вычисления производных. Целевая функция ожидаемо многоэкстремальна. Для повышения надежности результата в процесс встроена рандомизация.

2. 27-ТОЧЕЧНЫЙ ШАБЛОН ДЛЯ УРАВНЕНИЯ ГЕЛЬМГОЛЬЦА

Рассмотрим трехмерное уравнение Гельмгольца в однородной среде

Здесь $u(x,y,z)$ – неизвестная функция, $\Delta = \tfrac{{{{\partial }^{2}}}}{{\partial {{x}^{2}}}} + \tfrac{{{{\partial }^{2}}}}{{\partial {{y}^{2}}}} + \tfrac{{{{\partial }^{2}}}}{{\partial {{z}^{2}}}}$ – оператор Лапласа, $\omega $ – угловая частота, $c$ – (постоянная) скорость звука, $f(x)$ – функция источника (она исключается из дальнейших рассмотрений).

Введем в рассмотрение эквидистантную разностную сетку с шагами ${{h}_{1}},\;{{h}_{2}},\;{{h}_{3}}$. Узлами этой сетки являются точки с координатами $({{x}_{i}},{{y}_{j}},{{z}_{k}})$, пропорциональными соответствующим шагам сетки ${{x}_{i}} = i{{h}_{1}},$ ${{y}_{j}} = j{{h}_{2}},$ ${{z}_{k}} = k{{h}_{3}}$ с целыми коэффициентами $i,j,k$; значения $u({{x}_{i}},{{y}_{j}},{{z}_{k}})$ неизвестной функции в узлах сетки обозначаются ${{u}_{{ijk}}}$. Для конечно-разностной аппроксимации оператора $L$ уравнения (2.1) мы используем оператор ${{L}^{{(h)}}} = {{\Delta }^{{(h)}}} + \tfrac{{{{\omega }^{2}}}}{{{{c}^{2}}}}{{M}^{{(h)}}}$, где

Здесь ${{\Delta }^{{(1)}}}$, ${{\Delta }^{{(2)}}}$, ${{\Delta }^{{(3)}}}$ представляют собой аппроксимации оператора Лапласа разностными операторами с шаблонами, проиллюстрированными на фиг. 1. Там изображены три варианта шаблона для аппроксимации со вторым порядком второй производной ${{u}_{{zz}}}$ в центральной точке $(i,j,k)$ параллелепипеда: (1) иллюстрирует классическое трехточечное разностное отношение; на шаблоне (2) среднее арифметическое трехточечных отношений вдоль центральных линий граней; шаблон (3) использует ту же самую идею, но трехточечные разности берутся вдоль ребер. Скалярные весовые множители ${{\gamma }_{1}},{{\gamma }_{2}},{{\gamma }_{3}}$ подчинены условию ${{\gamma }_{1}} + {{\gamma }_{2}} + {{\gamma }_{3}} = 1$. Отметим, что положительность весов формально не требуется. Три варианта шаблонов для аппроксимации частной производной ${{u}_{{xx}}}$ и три варианта для ${{u}_{{yy}}}$ получаются аналогично.

Операторы ${{M}^{{(j,h)}}}$ представляют собой взятие среднего арифметического значений сеточной функции в точках, представленных черными шарами на соответствующих рисунках фиг. 2. Таким образом, в отличие от фиг. 1, шаблоны фиг. 2 полные. Подобно множителям ${{\gamma }_{j}}$ веса ${{w}_{j}}$ подчинены условию ${{w}_{1}} + {{w}_{2}} + {{w}_{3}} + {{w}_{4}} = 1$ и также не обязаны быть положительными.

Все рассмотренные шаблоны на фиг. 1 и 2 не выходят за пределы 27-точечного параллелепипеда. Соответственно шаблон оператора ${{L}^{{(h)}}}$ также получается 27-точечным.

3. ДИСПЕРСИОННАЯ ОШИБКА

Для дальнейшего нам понадобятся выражения для символов операторов, определяемых шаблонами фиг. 1 и фиг. 2. Для их получения мы действуем соответствующим разностным оператором на сеточную плоскую волну

Для шаблонов, изображенных на фиг. 1, получим В этих формулах $\xi = \mathop {\left( {{{\xi }_{1}},{{\xi }_{2}},{{\xi }_{3}}} \right)}\nolimits^t $, ${\text{|}}\xi {\text{|}} = \sqrt {\xi _{1}^{2} + \xi _{2}^{2} + \xi _{3}^{2}} $, символ ${{S}_{{{{\Delta }^{{(j,h)}}}}}}(\xi ;p)$ оператора ${{\Delta }^{{(j,h)}}}$, $j = 1,2,3$, связан с ним следующим соотношением В обозначениях символов явно прописана зависимость от параметра $p$, смысл которого объяснен ниже. Кроме того, использованы обозначения С учетом тождества $\tfrac{\omega }{{c{\text{|}}\xi {\text{|}}}} = 1$ выражения ${\text{|}}\xi {\text{|}}h$ приведены к виду ${\text{|}}\xi {\text{|}}h = 2\pi p$, где $p = 1{\text{/}}G$, а $G$ – это число точек на длину волны вдоль направления, определяемого наибольшим шагом сетки. Символ оператора ${{\Delta }^{{(h)}}} = {{\gamma }_{1}}{{\Delta }^{{(1,h)}}} + {{\gamma }_{2}}{{\Delta }^{{(2,h)}}} + {{\gamma }_{3}}{{\Delta }^{{(3,h)}}}$ есть линейная комбинация символов (3.1). Здесь $\gamma $ обозначает векторный параметр с компонентами ${{\gamma }_{1}},\;{{\gamma }_{2}},\;{{\gamma }_{3}}$.

Символ оператора ${{M}^{{(h)}}}$ является линейной комбинацией

символов В выражении ${{S}_{{{{M}^{{(h)}}}}}}(\xi ;p;w)$ (см. (3.2)) параметр $w$ обозначает вектор с компонентами ${{w}_{1}}$, ${{w}_{2}}$, ${{w}_{3}}$, ${{w}_{4}}$.

Дисперсионное соотношение для оператора ${{L}^{{(h)}}}$ принимает вид

Если параметры $\omega ,\;c,\;\xi $ конечно-разностной плоской волны подчинены соотношению (3.3), то отношение $\tfrac{\omega }{{{\text{|}}\xi {\text{|}}}}$ имеет смысл скорости ${{c}_{{num}}}$ распространения конечно-разностной волны. Соответственно выражение $\tfrac{\omega }{{c{\text{|}}\xi {\text{|}}}}$ представляет собой дисперсию скорости конечно-разностной волны, а разность $\tfrac{\omega }{{c{\text{|}}\xi {\text{|}}}} - 1$ можно назвать дисперсионной ошибкой. Для плоской волны, распространяющейся в направлении вектора $\xi $, формула для дисперсионной ошибки принимает вид

Формула (3.4) дает выражение для дисперсионной ошибки плоской волны, распространяющейся в направлении вектора $\xi $ и для выбранного значения параметра $p$. Используем равномерную метрику по параметрам $p$, $\xi $ для того, чтобы получить функцию

Оптимальные значения весов отыскиваются как решение задачи

В [5] вместо (3.5) используется функция

Определенным преимуществом такого подхода является гладкость функции $F$ и, как следствие, возможность использования классических методов оптимизации. В частности, градиент $F$ можно вычислить, дифференцируя интеграл по параметрам. Мы предпочитаем иметь дело с функцией $J$.

4. ВЫЧИСЛЕНИЕ ФУНКЦИИ $J(\gamma ,w)$

Графики дисперсионной ошибки, как функции сферических углов $\phi ,\;\theta $, параметризующих сферу ${\text{|}}\xi {\text{|}} = 1$, изображены на фиг. 3 для двух наборов весовых параметров. Видно, что дисперсионная ошибка кусочно-гладкая, имеет множество локальных максимумов и при значениях весов, близких к оптимальным, гладкость функции ухудшается.

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

Чтобы избежать громоздких выкладок по вычислению частных производных функции (3.4) аналитически, мы считаем их численно. Для этого организуем процесс вычисления последовательности разностных отношений на уменьшающихся приращениях аргументов до стабилизации значений этих отношений. Очевидной платой за это является несущественное увеличение времени счета.

При применении градиентных методов, имея направление градиента, нужно еще определить величину шага. На текущем значении шага мы вычисляем следующую точку и значение целевой функции в ней. Если значение целевой функции оказалось меньше, чем в текущей точке, новая точка не принимается, длина шага делится пополам и процесс повторяется с уменьшающимися шагами до тех пор, пока не получится найти большее значение или шаг не станет совсем маленьким, и будет найдена стационарная точка. Если же значение оказалось больше, то точка принимается, текущий шаг удваивается, в новой точке вычисляется градиент и процесс повторяется. Удвоение шага предусмотрено для того, чтобы в результате каких-то не очень удачных шагов не случилось так, что процесс идет с очень мелким шагом. Для того, чтобы не выйти за границы области по $p$, в случае если $p$-компонента вектора перемещения велика, мы укорачиваем длину шага так, чтобы не слишком сильно приблизиться к границе. Следует также учесть периодичность функции дисперсионной ошибки от сферических углов, добиваясь того, что точка всегда находится на основной карте $0 \leqslant \theta \leqslant \pi ,$ $0 \leqslant \phi < 2\pi $.

Градиентный метод оптимизации останавливается в точках локальных максимумов. Для поиска глобального максимума, который требуется в (3.5), мы применяем рандомизацию. Мы $K$ раз стартуем градиентный метод, используя в качестве стартовых точек случайно выбранные точки. Для получения надежных результатов значение $K$ должно быть достаточно большим. Экспериментально выяснено, что $K = 20$ – недостаточно, $K = 80$ – уже достаточно. При малых значениях $K$ решение задачи минимизации по весам $\gamma ,w$ получается нестабильным.

5. СИМПЛЕКС-МЕТОД

Для удобства читателя в этом разделе мы кратко излагаем элементы алгоритма безградиентной минимизации из [12], впервые предложенного в [11].

Рассматривается задача минимизации без ограничений вещественно-значной непрерывной функции $f({\mathbf{x}}),\;{\mathbf{x}} \in {{R}^{n}}$. Решение задачи осуществляется путем построения невырожденных симплексов, сходящихся к точке локального минимума. На $k$-м шаге симплекс ${{\Delta }_{k}}$ задан своим $n + 1$ вершинами ${\mathbf{x}}_{1}^{{(k)}},{\mathbf{x}}_{2}^{{(k)}}, \ldots ,{\mathbf{x}}_{{n + 1}}^{{(k)}}$.

Каждый шаг начинается с упорядочивания вершин ${{\Delta }_{k}}$ так, чтобы

Здесь $f_{j}^{{(k)}},\;1 \leqslant j \leqslant n + 1,$ обозначает $f({\mathbf{x}}_{j}^{{(k)}})$. Вершина ${\mathbf{x}}_{1}^{{(k)}}$ называется лучшей вершиной, вершина ${\mathbf{x}}_{{n + 1}}^{{(k)}}$ – худшей вершиной, вершина ${\mathbf{x}}_{n}^{{(k)}}$ – следующей за худшей. Подобным образом, значение $f_{{n + 1}}^{{(k)}}$ называется худшим значением функции, и т.п. Результатом $k$-го шага будет симплекс ${{\Delta }_{{k + 1}}}$, который соотносится с симплексом ${{\Delta }_{k}}$ одним из двух способов:

• ${{\Delta }_{{k + 1}}}$ и ${{\Delta }_{k}}$ имеют общую грань, образованную $n$ общими вершинами ${\mathbf{x}}_{1}^{{(k)}},{\mathbf{x}}_{2}^{k}, \ldots ,{\mathbf{x}}_{n}^{{(k)}}$;

• ${{\Delta }_{{k + 1}}}$ и ${{\Delta }_{k}}$ имеют одну общую вершину ${\mathbf{x}}_{1}^{{(k)}}$ (лучшая вершина в ${{\Delta }_{k}}$).

Оригинальный критерий остановки процесса, предлагаемый авторами, выглядит так

где ${\mathbf{\bar {x}}} = \sum\nolimits_{i = 1}^n \,{{{\mathbf{x}}}_{i}}{\text{/}}n$ – центр $(n - 1)$-мерной грани симплекса. Мы дополняем его еще тремя условиями – в нашем варианте вычисления продолжаются, пока выполняются следующие четыре условия.

1. Среднеквадратичное уклонение значений функции в точках симплекса превышает установленный порог (см. (5.2)).

2. Отношение объема текущего симплекса к объему исходного симплекса превышает установленный порог.

3. Обусловленность матрицы, определяемой ребрами симплекса, выходящими из худшей вершины, не превышает установленный порог.

4. Число сделанных шагов не превосходит некоторого предела.

Точки, в которых симплекс-метод останавливается по нарушению первого условия, претендуют на то, чтобы быть стационарными точками функции $f({\mathbf{x}})$ (локальными минимумами). Для того, чтобы приблизиться к глобальному минимуму в процесс минимизации, мы включаем рандомизацию. Вновь найденная стационарная точка получает небольшое случайное возмущение, и процесс вновь стартует с новой начальной точки. После следующей остановки нужно сравнить значение целевой функции в новой стационарной точке с предыдущим значением.

Несмотря на огромную популярность симплекс-метода, его сходимость может вызывать нарекания (см., например, [13]).

6. ЧИСЛЕННЫЕ ПРИМЕРЫ

В нашем первом примере все шаги сетки одинаковы, интервал по $p$ от 0.01 до 0.2, что соответствует условию “пять и более точек на длину волны”. Стартовое значение весовых параметров $\gamma = (0.6401,0.3285,0.0314)$, $w = (0.5365,0.3276,0.2051, - 0.0692)$, дисперсионная ошибка в стартовой точке равна $1.67 \times {{10}^{{ - 3}}}$.

В процесс вычисления $J(\gamma ,w)$ мы включали рандомизацию, 80 раз ($K = 80$) стартуя градиентную максимизацию (по параметру $p$ и направлениям распространения плоской волны) дисперсионной ошибки. Включая рандомизацию в процесс минимизации $J(\gamma ,w)$, мы пять раз запускали симплекс-метод для того, чтобы минимизировать функцию $J(\gamma ,w)$. Как это описано в конце раздела, в качестве новой стартовой точки выбиралось некоторое случайное возмущение полученной стационарной точки. Результаты оптимизации представлены в пяти столбцах табл. 1. Значения в этой таблице приведены только с 8 значащими цифрами, но заметим, что в реальности весовые коэффициенты из разных столбцов совпадают с точностью до 16 знаков, значения целевой функции – с 10 знаками. Число шагов симплекс метода от стартовой до стационарной точки приведено в строке ${{n}_{{iter}}}$. В последней строке таблицы расположен номер критерия, по которому был остановлен симплекс метод. Во всех столбцах там располагается число 1, иначе говоря, процесс остановился в стационарной точке.

На фиг. 9 изображены 100 кривых дисперсионной ошибки в зависимости от параметра $p$, соответствующих 100 случайным векторам распространения плоских волн. Ошибка вычислялась для оптимизированной разностной схемы с весами из табл. 1. Легко видеть, что для разностных сеток с числом точек на длину волны, равным пяти и больше, дисперсионная ошибка не превосходит 0.1%.

Наш второй численный пример касается разностной сетки с отношением шагов $1:3:5$. мы стартовали оптимизацию с классической семиточечной схемы $\gamma = (1,0,0),w = (1,0,0,0)$, используя $K = 100$. Как и раньше, полученная точка выбиралась в качестве стартовой для следующего шага. Можно вновь заметить стабилизацию результатов.

Кривые дисперсионной ошибки в зависимости от параметра $p$ для оптимизированной разностной схемы с весами из пятого столбца табл. 2 представлены на фиг. 10. Кривые построены для 100 случайно сгенерированных направлений распространения плоской волны. Как следует из фиг. 10, для разностных сеток с числом точек на длину волны, равным пяти и больше, дисперсионная ошибка не превосходит 0.1%.

Кроме численной дисперсии, в решениях разностного уравнения также присутствует явление анизотропии, проиллюстрированное на фиг. 11. Дисперсионная ошибка оптимизированной схемы представлена в виде трехмерной поверхности, гомеоморфной сфере. Для ее получения соответствующий радиус-вектор был умножен на множитель $1 + a \cdot d$, где $d$ – дисперсионная ошибка в заданном направлении, $a$ – коэффициент растяжения. Для сравнения приведена также единичная сфера. Точки пересечения поверхности со сферой определяют направления с нулевой дисперсионной ошибкой, точки снаружи сферы соответствуют направлениям, для которых скорость распространения разностной плоской волны превышает истинную скорость звука. Соответственно точки внутри сферы относятся к направлениям, вдоль которых плоские разностные волны распространяются со скоростью, меньшей чем истинная скорость звука.