Акустический журнал, 2021, T. 67, № 5, стр. 533-541

ВВЕДЕНИЕ

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

В настоящее время в литературе представлены несколько групп методов расчета трехмерных звуковых полей, среди которых можно отметить геометроакустические (лучевые) методы [1, 2], теорию гауссовых пучков [3] и метод параболического уравнения [4, 5]. Прикладные программы, основанные на методах из первых двух групп, представлены в сети Интернет [6], однако их применение в случае трехмерных волноводов мелкого моря с неоднородной батиметрией и гидрологией затруднено нетривиальным поведением семейства выходящих из источника лучевых траекторий [1].

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

Настоящая работа представляет собой описание комплекса программ AMPLE, предназначенного для решения задач расчета акустических полей в волноводах мелкого моря весьма общего вида [7], разработанного в ТОИ ДВО РАН и основанного на численном решении широкоугольных модовых параболических уравнений (ШМПУ) [8, 9]. Под задачей общего вида далее понимается задача, в которой рельеф дна задается произвольной кусочно-гладкой функцией общего вида $h = h\left( {x,y} \right)$ (например, заданной своими значениями на некоторой сетке $\left( {{{x}_{i}},{{y}_{j}}} \right)$), поле скорости звука определяется набором вертикальных разрезов (ВРСЗ) в нескольких точках акватории $c = c\left( {z,{{x}_{s}},{{y}_{s}}} \right)$ (все они считаются кусочно-гладкими функциями глубины z, например, тоже заданными на некоторой сетке), а параметры дна (скорость звука ${{c}_{b}}$, плотность ${{{{\rho }}}_{b}}$ и поглощение ${{{{\beta }}}_{b}}$) имеют произвольные значения либо произвольным образом зависят от глубины z. Комплекс программ AMPLE [7] при этом позволяет выполнить расчет трехмерного поля акустического давления $p\left( {x,y,z} \right)$, формируемого точечным источником тонального акустического сигнала заданной частоты f. Расчет поля выполняется на заданной трехмерной сетке $\left( {{{x}_{n}},{{y}_{m}},{{z}_{k}}} \right)$. Мы полагаем, что любой комплекс программ, ориентированный на использование в практических задачах, должен обеспечивать данный функциональный минимум, который, разумеется, может быть и значительно расширен.

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

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

ШИРОКОУГОЛЬНЫЕ МОДОВЫЕ ПАРАБОЛИЧЕСКИЕ УРАВНЕНИЯ И ИХ ЧИСЛЕННОЕ РЕШЕНИЕ

В данном разделе мы приводим краткое описание методики решения модовых параболических уравнений, использованной в нашем комплексе программ. Заметим, что на настоящий момент реализовано только дробно-линейное приближение (аппроксимация Паде 1-1) для операторного квадратного корня. В настоящее время завершается разработка новой версии, в котором для решения ШМПУ будет использоваться метод SSP [10, 11] с произвольным порядком аппроксимации Паде для экспоненты операторного квадратного корня. Как следует из формул, приведенных в [9], выражения, по которым в маршевой схеме осуществляется пересчет от x_n к x_{n + 1}, при этом практически не изменяются. Подробное описание всевозможных методов численного решения ШМПУ представлено в работе [9], и краткое описание ниже приведено лишь для формирования у читателя полной картины без необходимости обращаться к публикациям математического характера.

Поле акустического давления $p\left( {x,y,z} \right)$ (где $z$ обозначает глубину, а $x,y$ – координаты горизонтальной плоскости), создаваемое точечным источником в трехмерном волноводе мелкого моря, расположенным в точке $x = y = 0,z = {{z}_{s}}$ и имеющим частоту $f$, описывается трехмерным уравнением Гельмгольца [10, 12]

где ${{\rho }}\left( z \right)$ – плотность среды, $c\left( z \right)$ – скорость звука, ${{\omega }} = 2{{\pi }}f$ – циклическая частота, ${{\beta }}$ – коэффициент поглощения, ${{\eta }} = \frac{1}{{40{{\pi }}{{{\log }}_{{10}}}e}}$. Решение такого уравнения в адиабатическом приближении может быть представлено в виде модового разложения [10, 12] где ${{A}_{j}}\left( {x,y} \right)$ – модовые амплитуды, удовлетворяющие уравнениям горизонтальной рефракции [8, 10] ${{{{\varphi }}}_{j}}\left( {z,x,y} \right)$ – модовые функции, а ${{k}_{j}}\left( {x,y} \right)$ – соответствующие им волновые числа, которые могут быть получены из решения акустической спектральной задачи [10].

Уравнение (3) может быть преобразовано в широкоугольное модовое параболическое уравнение

где ${{k}_{{j,0}}}$ – опорное волновое значение, ${{{{\alpha }}}_{0}} = a - 1,$ ${{{{\alpha }}}_{1}} = b - c$, $a,b,c$ – коэффициенты дробно-линейной аппроксимации оператора квадратного корня [10]

Решение уравнения (4) производится с использованием численной схемы, основанной на стандартном методе Крэнка–Николсон, на равномерной прямоугольной сетке

Результат решения при этом представляет собой двумерный массив значений $\mathcal{A}_{j}^{{n,m}} \approx {{\mathcal{A}}_{j}}\left( {{{x}_{n}},{{y}_{m}}} \right).$

Решение ШМПУ всегда рассматривается в неограниченной в направлении y области, что приводит к необ ходимости искусственного ограничения расчетной области. Чтобы в таком случае эмулировать распространение звука в неограниченной области, были использованы дискретные граничные условия прозрачности, изначально разработанные для “вертикальных” параболических уравнений [13]. Условия устанавливаются на левой и правой границах области при $y = {{y}_{0}}$ и $y = {{y}_{M}}$. При использовании таких граничных условий необходимо дополнительное упрощающее предположение о том, что волновые числа ${{k}_{j}}$ на обеих границах области не зависят от продольной (эволюционной) переменной x, то есть

В завершение данного раздела еще раз отметим, что все расчеты в рассматриваемом комплексе программ выполняются в рамках адиабатического приближения. Кроме того, переход от уравнения горизонтальной рефракции (3) к ШМПУ (4) предполагает, что мы пренебрегаем обратным рассеянием звука в направлении x. Наконец, дробно-линейная аппроксимация Паде характеризуется эффективной апертурой около ±15° в горизонтальной плоскости.

Часть этих ограничений будет снята при дальнейшем совершенствовании комплекса программ. Так, замена (4) на псевдодифференциальное МПУ (ПДМПУ) [8, 9] позволяет практически устранить последнее ограничение и добиться апертуры в горизонтальной плоскости, составляющей по меньше мере ±75°. Однако, по нашему опыту, дробно-линейная аппроксимация вполне достаточна для решения почти любых практических задач (достаточно сложно предложить пример задачи, в которой в поле на заданной акустической трассе существенный вклад будут вносить лучи, выходящие из источника под углами более 15° к этой трассе). Использование адиабатического приближения в широком классе практических задач акустики мелкого моря оправдано тем, что на шельфе типичные углы наклона дна обычно не превышают 1°. Тем не менее, при моделировании распространения звука с шельфа в глубокий океан взаимодействие мод начинает играть решающую роль, и потому в настоящее время нами разрабатывается обобщение теории ПДМПУ, которое позволит учесть данный эффект.

РАСЧЕТ АКУСТИЧЕСКИХ МОД

Первым шагом при расчете акустического поля по формуле (2) является вычисление модовых функций ${{{{\varphi }}}_{j}}\left( {z,x,y} \right)$ и волновых чисел ${{k}_{j}}\left( {x,y} \right)$ на некоторой сетке $\left( {{{x}_{n}},{{y}_{m}},{{z}_{k}}} \right)$. Это вычисление осуществляется путем многократного численного решения акустической спектральной задачи [10].

Для расчета акустических мод на языке C++ был разработан программный комплекс CAMBALA (сокр. от Coupled Acoustic Modes with Bottom Attenuation in Linear Acoustics) [14]. Реализованы два варианта работы CAMBALA. В первом варианте это отдельная программа, которая получает на вход файл сценария, проводит необходимые вычисления и записывает результаты в выходные файлы (описание ввода, вывода и использования программы CAMBALA в этом варианте приведено в соответствующем руководстве пользователя) [14]. Во втором варианте это библиотека, функции которой могут быть вызваны из внешней программы. Именно этот вариант используется нами как часть комплекса программ для расчета акустического поля по формуле (2).

Основой CAMBALA является класс NormalModes. Для решения спектральной задачи в данной точке волновода программа, использующая CAMBALA в качестве библиотеки, должна сформировать объект данного класса и передать информацию о зависимости параметров волновода от глубины z в данной точке $\left( {{{x}_{n}},{{y}_{m}}} \right)$ путем присвоения значений соответствующим полям класса, указанным в заголовочном файле. Зависимости параметров волновода от глубины считаются кусочно-линейными функциями, определяемыми набором слоев, параметры которых заданы в четырех массивах равной длины, совпадающей с числом этих слоев (каждый из них представляет собой поле в классе NormalModes). Кроме этого, NormalModes имеет поля, с помощью которых задаются горизонты для расчета значения модовых функций, частота звука, количество мод, количество точек на 1 м, используемое при дискретизации спектральной задачи, порядок экстраполяции Ричардсона и подмножество мод, которые необходимо рассчитать. Имеется также возможность указать алгоритм решения спектральной задачи. В настоящий момент возможно использовать два алгоритма. В обоих случаях выполняется дискретизация акустической спектральной задачи, которая приводит ее к задаче на собственные значения некоторой трехдиагональной матрицы [10, 15]. По умолчанию используется алгоритм “alglib”, в рамках которого для вычисления собственных значений этой матрицы будет использоваться библиотека ALGLIB [16], в которой реализован метод Штурма (с хорошим уровнем оптимизации вычислений). Альтернативой ему является алгоритм “spectra”, реализованный в одноименной библиотеке [17], основанной на методе итераций Арнольди.

Исходный код CAMBALA является платформонезависимым, причем при его разработке использован стандарт языка C++17 [14] (таким образом, в настоящее время CAMBALA, по-видимому, является единственной библиотекой для расчета акустических мод, разработанной в рамках современных стандартов программирования). Для сборки доступны два проекта Visual Studio, первый из которых позволяет собрать как 32-битный, так и 64-битный exe-файл для операционной системы семейства Windows. Второй проект Visual Studio предназначен для сборки библиотеки. Набор инструкций файла Makefile позволяет собрать 32-битный либо 64-битный исполняемый файл (либо соответствующую библиотеку) в Linux-подобных операционных системах. Для удобства сборки также имеются файлы конфигурации в формате CMake, автоматически загружающие все необходимые зависимости.

Разработанный программный комплекс CAMBALA размещен в свободном доступе в сети Интернет [14]. В данном github-репозитории находятся: исходный код, 64-битные исполняемые файлы для Windows и Linux, 64-битные библиотеки для Windows и Linux, руководство по сборке и запуску, примеры сценариев, а также примеры запуска и подробная документация.

АЛГОРИТМ РАСЧЕТА ПОЛЯ

Расчет акустического поля по формуле (2) производится программой AMPLE [7]. Параметры вычислений задаются в конфигурационном файле в формате JSON, представляющим собой объект, содержащий список пар ключ/значение. Конфигурационный файл позволяет задавать параметры среды: границы и количество точек вычислительной области, описание параметров слоев воды и дна, данные гидрологии и батиметрии, частоты источника, тип и параметры начальных условий, скалярные параметры класса NormalModes, и др. В репозитории также содержится более подробное описание параметров и примеры конфигурационных файлов.

Данная программа инициирует вычисление ${{k}_{{j,0}}}$ и ${{{{\varphi }}}_{j}}\left( {{{z}_{s}},0,0} \right)$, используя библиотеку CAMBALA и создавая объект класса NormalModes. После этого строятся начальные условия ${{\mathcal{A}}_{j}}\left( {0,y} \right)$ задачи Коши для ШМПУ (4). Далее производится расчет модовых функций ${{{{\varphi }}}_{j}}\left( {x,y} \right)$ и соответствующих им волновых чисел ${{k}_{j}}\left( {x,y} \right)$. Расчет производится с использованием пакета CAMBALA на основе параметров среды: батиметрии, гидрологии, структуры дна и т.д. Этот расчет проводится на сетке, отличной от вычислительной и зачастую более грубой, для борьбы с этим используется линейная интерполяция. Следующим выполняется вычисление решения ШМПУ с граничными условиями прозрачности. В дальнейшем вычисление строится следующим образом:

• производится запуск $w + 1$ потока выполнения;

• первые $w$ потоков выполняют вычисление ${{\mathcal{A}}_{j}}$ для соответствующих интервалов модовых функций и волновых чисел $j_{w}^{0} \leqslant j \leqslant j_{w}^{1}$, используя метод прогонки, при этом хранится только предыдущая строка поля $y = {{y}_{{m - 1}}}$, а каждая новая строка $y = {{y}_{m}}$ передается в специальную callback функцию для последующей обработки;

• $w + 1$ поток выполнения ожидает рассчитанные остальными потоками строки решения, производит суммирование по формуле $\left( 2 \right)$ и передает результат в callback функцию вывода решения ШМПУ.

• Вычислительная сложность может быть оценена как $O\left( {F + JN\left( {N + M} \right)} \right)$, где $F$ – время, затраченное на вычисление модовых функций и волновых чисел. Вычисление одной строки решения требует $O\left( N \right)$ для вычисления граничных условий $\left( {12} \right)$ и $O\left( M \right)$ для решения системы методом прогонки. Такие вычисления повторяются $N$ раз для $J$ мод. Так как в процессе вычислений хранится только предыдущая строка решения и значения функции ${{\mathcal{A}}_{j}}$ на границах вычислительной области, затраты памяти составляют $O\left( {J\left( {N + M} \right)} \right)$.

Результаты вычислений сохраняются в папке, указанной в параметрах командной строки. Папка содержит содержит следующие файлы:

• config.json – описание всех параметров, использованных во время вычислений;

• meta.json – описание выходных файлов и папок, параметров командной строки, время, затраченное на проведение вычислений.

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

ТЕСТОВЫЕ ПРИМЕРЫ

Репозиторий комплекса программ включает в себя набор тестовых примеров, которые одновременно могут быть использованы в качестве образца при его конфигурировании на практике. В первом и втором примерах решаются задачи распространения звука в идеализированных модельных волноводах. Тестовый пример 1 представляет собой классическую задачу о распространении звука в прибрежном клине с наклоном 2.86° и проницаемым дном [18] (см. рис. 1), в которой решение может быть найдено аналитически с использованием метода изображений. Мы не будем перечислять все параметры волновода в данном примере, поскольку он используется для тестирования различных методов моделирования звуковых полей уже более 20 лет [3–5, 9, 10]. Заметим лишь, что они указаны в соответствующем конфигурационном файле, который имеется в репозитории AMPLE [7]. В данном примере принято рассчитывать уровни акустического поля (формируемого точечным источником с частотой 25 Гц) как функцию горизонтальных координат (x, y) на горизонте z = 30 м. Это поле, рассчитанное с помощью комплекса AMPLE, представлено на рис. 2. На странице с описанием комплекса [7] имеется скрипт MATLAB, который открывает файлы вывода AMPLE, читает оттуда рассчитанные поля и строит по ним рисунки, подобные представленным в данном разделе. При построении рисунков также выполняется сопоставление результата, полученного с помощью AMPLE, с результатом, рассчитанным с использованием метода изображений. При этом наблюдается практически идеальное совпадение графиков уровня поля. Мы полагаем, что авторам всех методов расчета звуковых полей и основанных на них комлексов программ необходимо сопровождать их данным примером. Отметим, что в данном примере поле, рассчитанное в рамках приближений метода адиабатических ШМПУ, не отличимо от поля, рассчитанного с помощью метода изображений [9].

Во втором тестовом примере рассматривается распространение звука в волноводе мелкого моря с подводным каньоном (см. схему на рис. 3), наличие которого обуславливает фокусировку звука в горизонтальной плоскости в окрестности его оси [9, 19]. Эталонное решение для примера 2 также представлено на странице AMPLE. Оно рассчитано с использованием метода виртуальных источников [20]. Рельеф дна в данном примере задается функцией

где h₀ = 20 м, Δh = 15 м, σ = 7 × 10⁻⁴ м⁻¹. Дно считается жидкой средой со скоростью звука c_b = = 1800 м/с, плотностью ρ_b = 2 г/см³ и поглощением β = 0.5 дБ/λ. Мы рассчитываем поле, формируемое точечным источником с частотой 150 Гц, расположенным в точке x = 0, y = 0, z = 10 м. Результатом решения задачи является горизонтальный разрез этого поля плоскостью z = 10 м (см. рис. 4). Как показано в [9], несмотря на перепады глубины, обеспечивающие отсечку двух из четырех водных мод, формируемых источником звука, решение, полученное с помощью адиабатических ШМПУ, практически идеально совпадает с эталонным. Мы полагаем, что данный пример также может стать обязательным элементом тестирования любых методов моделирования звуковых полей в мелком море.

Третий тестовый пример, конфигурационный файл для которого имеется в репозитории AMPLE [7], представляет собой практическую задачу распространения звука на акватории шельфа о. Сахалин. Батиметрические данные для данной акватории были взяты нами из базы данных SRTM (см. рис. 5), а зависимость скорости звука от глубины получена путем аппроксимации данных натурных измерений, выполненных на этой акватории в летний период с помощью сигмоиды (см. рис. 6). Разумеется, для такой задачи эталонного решения не существует. Тем не менее, с нашей точки зрения, пользователям программы будет удобно иметь конфигурацию для запуска ее решения, включающую в себя близкие к натурным данные о гидрологии и батиметрии. В данном примере также рассчитывается горизонтальный разрез звукового поля (на глубине z = 10 м), формируемого точечным источником с частотой 200 Гц. Результаты расчета представлены на рис. 7.

Еще раз отметим здесь, что первые два примера в данном случае служат для апробации не математического метода ШМПУ (которая выполнена нами ранее [9]), а конкретной его программной реализации.

ЗАКЛЮЧЕНИЕ

В настоящей работе представлено описание комплекса программ AMPLE, предназначенного для численного моделирования акустических полей в мелком море и основанного на методе адиабатических ШМПУ. Решение тестовых задач показывает, что, несмотря на ограничения, налагаемые на точность расчетов адиабатическим приближением, а также приближением однонаправленного распространения, данный комплекс и лежащая в его основе методика позволяют выполнять достаточно точные для многих приложений расчеты звуковых полей в мелком море с неоднородным рельефом дна.

По-видимому, AMPLE является первым открытым отечественным программным продуктом, предназначенным для моделирования трехмерных звуковых полей. При его разработке использованы современные стандарты программирования, что отличает его от известных комплексов KRAKEN, ORCA, BELLHOP, авторы которых в основном придерживаются весьма консервативного стиля написания программ. Отметим также, что число общедоступных комплексов программ, обеспечивающих надежность результатов расчета трехмерных звуковых полей хотя бы на нескольких модельных примерах невелико. В этой связи мы надеемся, что AMPLE окажется востребованным специалистами из разных областей акустики океана. К числу приложений, в которых AMPLE уже используется в настоящее время в ТОИ ДВО РАН, относятся задачи определения зон акустической тени (в горизонтальной плоскости), которые могут возникать вследствие горизонтальной рефракции звука при размещении источников навигационных сигналов на шельфе (см., например, [21]), а также задачи прогнозирования уровней антропогенных шумов при проведении акустического мониторинга акваторий [22].

В дальнейшем предполагается развитие комплекса программ AMPLE при сохранении основных принципов его функционирования и способов представления входных и выходных данных. Среди ближайших планируемых усовершенствований – добавление возможности расчета модовых амплитуд с использованием метода SSP (split-step Padé, метод аппроксимации экспоненты операторного квадратного корня [9–11]). В настоящий момент также разрабатывается методика, которая позволит выполнять расчеты с учетом взаимодействия мод [10] (т.е. окончательно снять ограничения, налагаемые адиабатическим приближением). Предполагается также добавление возможности моделирования импульсных сигналов и расчета их интегральных энергетических характеристик.

Данная работа выполнена в рамках тем гос. задания ТОИ ДВО РАН “Моделирование разномасштабных динамических процессов в океане” (0211-2021-0009). Работа А.Т. и П.П. также частично поддержана Эксон Нефтегаз Лимитед.