Содержание к диссертации
Введение
1. Введение и обзор проблемы 4
1.1. Постановка задачи 6
1.2. Слоистый океан ... 7
1.3. Двумернонеоднородный океан 17
1.4. Трехмернонеоднородный океан 19
2. Слоистый океан 21
2.1. Граничные условия на дне океана 22
2.2. Непоредственная численная оценка интегрального представления поля 32
2.3. Метод нормальных волн 42
2.3.1. Вычисления собственных функций путем разложения в ряд... 42
2.3.2. Вычислние собственных функций комбинированным методом. 51
2.3.3. Некоторые численные примеры 58
2.3.4. Некоторые интегральные соотношения для нормальных волн 66
2.3.5. Модовая структура поля в слоистом океане 70
2.3.6. Резонансное взаимодействие нормальных волн с тонкой структурой скорости звука в слоистом океане 90
2.4. Элементы теории распространения звука в слоистом океане в терминах нового асимп тотического представления решений 95
2.4.1. Модифицированный метод ВКБ 96
2.4.2. Поперечная функция Грина 104
2.4.3. "Лучевое" представления поля 112
3. Двух- и трежернонеодеюродный океан 119
3.1. Постановка задачи 120
3.2. Аппроксимация формы дна океана 122
3.3. Аппроксимация свойств среды (первый способ)... 134
3.3.1. Коэффициент преломления 135
3.3.2. Лучевые уравнения 138
3.4. Аппроксимация свойств среды (второй способ).. 143 3.4.1. Барицентрические и "слоистые" координаты 144
3.5. Метод суммирования гауссовых пучков (МСГП)... 154
3.5.1. Асимптотика поля точечного гармонического источника 155
3.5.2. Аналитические конструкции гауссова пучка при некоторых аппроксимациях скорости звука ... 162
3.5.3. Численная реализация МСГП 165
3.6. Примеры расчетов поля лучевым методом и МСГП 168
4. Некоторые обратные задачи в слоистом океане 181
4.1. Асимптотическая связь спектра нормальных волн в слоистом океане с завис шлостью скорости звука от глубины 181
4.2. Определение акустической неоднородности среды с помощью звуковых сигналов 187
5. Основные результаты и заключение 198
6. Литература
- Двумернонеоднородный океан
- Вычисления собственных функций путем разложения в ряд...
- Аппроксимация свойств среды (первый способ)...
- Определение акустической неоднородности среды с помощью звуковых сигналов
Введение к работе
В подводной акустике за последние 15-20 лет сформировалось целое направление исследований - математическое моделирование процессов распространения звуковых сигналов в океане. Постановка экспериментов по распространению звука в океане требует больших затрат времени и денег, при этом практически невозможна постановка "чистых" экспериментов, в силу большого числа неконтролируемых факторов, таких как состояние среды, поверхности, структуры морского дна и так далее. Численные эксперименты неизмеримо дешевле, быстрее и возможны в широком диапазоне полностью контролиуемых условий распространения звука. С усовершенствованием программ вычисления полей и с повышением быстродействия ЭВМ становится реальным создание таких численных моделей, которые дают значения звуковых полей и элементов их структур быстрее, чем это происходит в реальных физических ситуациях, что может оказать решающее влияние на решение многих важных прикладных и исследовательских задач. Численные эксперименты дают возможность полностью проанализировать структуру звуковых полей в пространстве и времени как результат тех или иных условий распространения, местоположения и конфигурации источников и приемников звука, частот излучения и так далее.
Математическая формулировка задач о распространении звука в различных средах сложилась довольно давно, однако численная реализация решений уравнений движения возможна с помощью различных средств, идущих как от классической математики (асимптотика, интегральные преобразования и т.п.), так и от вычислительной математики (разностные схемы, сплайны и т.п.). По глубокому убеждению автора наиболее перспективными являются такие методы расчета звуковых полей (независимо от их происхождения),которые обладают с одной стороны высоким быстродействием и с другой сторо ны являются "открытыми", то есть позволяют ввести в задачу дополнительные усложнения, обусловленные учетом новых физических факторов. Очень важным свойством алгоритмов является также простота их структуры, позволяющая создавать комплексы программ, нацеленные на решение широкого круга прикладных задач. Быстродействие и "открытость" алгоритмов часто являются взаимоисключающими свойствами, так как быстродействие алгоритма почти всегда основано на использовании тончайших нюансов задачи, что приводит к построению логически сложных и трудно поддающихся перестройке программ. С другой стороны "открытость" или простота структуры алгоритма часто идет в ущерб быстродействию. Разумеется решающее влияние на выбор того или иного метода решения задачи и на структуру программ оказывают возможности вычислительных средств, в первую очередь длина мантиссы, наличие быстрой внешней памяти на дисках, наличие спецпроцессоров и так далее. Изложенные ниже аналитические и численные подходы к задачам распространения звука демонстрируют попытки автора учесть указанные факторы с целью создания эффективных методов вычисления звуковых полей.
Двумернонеоднородный океан
Если вычислять значение функции Грина при равноотстоящих значениях т таких, что то вычисление функций Бесселя, требующее больших затрат машинного времени, можно произвести лишь для начальных значений (Ук)0 , (Ук){ » а значение остальных получить с помощью рекуррентных соотношений для Бесселевых функций. Авторы провели дополнительное ис-ледование рекуррентных соотношений в области, где последние неустойчивы, то есть при ЇК І КІҐУІІ ЇК + І Поскольку величины Нк » определяющие шаг по , различны в различных слоях, шаг выбирается по формуле где Нтах наименьшее общее кратное всем Н к . Полученная таким образом функция Грина подставляется в интегральное представление (I) и последний интеграл оценивается с помощью алгоритма быстрого преобразования Фурье, после замены функции Ханкеля ее асимптотическим разложением.
Аналитическое совершенство алгоритма ррр имеет оборотную сто рону - большой объем памяти ЭВМ, возрастающий при прочих равных условиях с частотой излучаемого звука. В книге [65] указано, что для относительно небольших частот излучения в океане средней глубины (70 гц, 3 км) алгоритм требует более 50 к слов машинной памяти, что делает его недоступным для подавляющего большинства отечественных машин и микро-ЭВМ. В связи с этим нами был разработан альтернативный алгоритм быстрой оценки звукового поля в слоистом океане (см. Главу 2), требующий меньше памяти и обладающий достаточно высоким быстродействием.
Большой набор методов вычисления звукового ПОЛЯ в слоистом океане предназначен не только и не столько для вычисления самих звуковых полей, а для численного исследования разнообразных физических эффектов, наблюдаемых при распространении звука. В работах [6,6-77] рассматривалась пространственно-частотная интерференционная структура звукового поля в слоистом океане. В Главе 2 дан теоретический анализ этого явления и приведены численные расчеты, подтверждающие результаты теории. При этом существенно используется асимптотическая связь между лучами и нормальными волнами в слоистой среде, для которой используется соотношение, связывающее длину цикл луча, номер моды, фазовую и групповую скорости и частоту излучения [78].
В работах [79-81] обнаружен эффект смещения зон конвергенции в реальных океанических экспериментах по сравнению с лучевыми расчетами. В Главе 2 дано описание одного из возможных механизмов, вызывающих указанное смещение.
В случае, когда скорость звука в океане, либо форма морского дна зависят от продольной координаты, число эффективно работающих методов вычисления звукового поля резко уменьшается, по сравнению с рассмотренным выше случаем слоистого океана. Существуют единичные примеры точного решения задачи Ї1 и все методы расчета поля являются приближенными.
Метод нормальных волн, при переходе к неслоистой среде, преобразуется в так называемый метод поперечных сечений [1,82-86] . В этом методе решение исходного уравнения Гельмгольца отыскивается в виде разложения по локальным собственным функциям поперечного оператора с неизвестными коэффициентами. Подстановка этих разложений в уравнение Гельмгольца приводит к бесконечной системе обыкновенных дифференциальных уравнений для искомых коэффициентов. В адиабатическом приближении считается, что матрица этой системы диагональна и для амплитуд мод получаются обыкновенные дифференциальные уравнения.
Условия применимости адиабатического приближения являются очень жесткими [1,88] и поэтому были предприняты попытки учесть в той или иной степени взаимодействие между модами [89,90] . В [91-100] выписаны уравнения метода поперечных сечений для различных физических задач о распространении звука в неслоистых волноводах и приведены аналитические и численные результаты исследования их решений.
Формализм лучевого метода с легкостью переносится на двумерно-неоднородный океан, однако, только в случае специально подобранных способов аппроксимации коэффициента преломления среды удается получить куски лучей в аналитической форме, чтобы повысить быстродействие алгоритмов. В [Ю2] использована блочная структура аппроксимации коэффициента преломления с аналитической формой луча в каждом слое. Ццея Вагина [103,104] состояла в том, чтобы искать форму луча в виде разложения по степеням длины дуги, что привело к созданию высокопроизводительного алгоритма расчета как в слоистой, так и в неслоистой средах. В [105,106] получено дифференциальное уравнение для интенсивности в неслоистой среде. В [107] показано, что для того чтобы избежать ложных каустик для этого случая, необходима дважды непрерывно дифференцируемая аппроксимация коэффициента преломления. Большая группа работ [I09-II4] посвящена лучевым алгоритмам и расчетам по ним в двумернонеоднородном океане, в том числе и с неоднородностями, порождаемыми океанскими течениями и приливами. Интересный вариант лучевого метода, полученный с помощью разложения Гамильтониани задачи в окрестности луча в квадратичную форму изложен в [115,116] .
Вычисления собственных функций путем разложения в ряд...
Поле точечного гармонического источника звука в слоистом оке ане, то есть в цилиндрических коэффициентах в области 0; О Z /-/ -Л Lp лежащей на упругом полупространтсве со ско ростями Ст - Ст(г); Сб- Се (?) есть [і] где tfnfjn z) есть решение задачи на собственные значения и " + f«2 -тг7и? л (2) примем Q f? ) в зависимости от соотношения между фазовой скоростью гь-Си/ и скоростями Ст (н) и Сь(н) задается формулами (2.1.38), (2.1.42) и (2.1.46). В последующих разделах мы остановимся на различных способах вычисления собственных функций Уп. , включающих в себя как традиционные подходы (метод ВКБ, разложение в ряд по степеням ii ), так и разработанные автором модифицированный метод ВКБ.
Вычисление собственных функций путем разложения в ряд. Метод вычисления решений уравнения ( 2-) путем разложения в ряд по глубине Z впервые был применен в работах / 7,8j. Однако в этих работах не была исследована ошибка, возникающая при суммировании ряда на конкретной ЭВМ и не был построен устойчивый метод интегрирования. В настоящем разделе изложен вариант метода, восполняющий указанные недостатки.
Поскольку в эксперименте зависимость скорости звука от глу - 43 бины 0 (г) обычно задается в виде таблицы значений CL(2i), і-1,...ґі9 то между этими значениями либо скорость, либо коэффициент преломления обычно аппроксимируется полиномом. В результате мы получаем на каждом отрезке / ; f i_hi 7 вместо уравнения ( %, ) уравнение с полиномиальными коэффициентами и, естественно, для его интегрирования воспользоваться полиномиальным разложением решения на каждом отрезке [z ; сн1 Здесь приведен пример простейшей аппроксимации, когда при Zefe yj
Последнее, очевидно, сводится к уравнению Зйри вдали от точки поворота, однако мы будем производить прямое интегрирование этого уравнения по соображениям, которые будут выяснены к концу настоящего раздела.
Зададим на поверхности 2 О условия Пі , 4(ю) = (4) и на дне, при - Ц » условия r(iu)=-6Q(-i)} Ъ(1н)=ё (5) где CL и 8 - произвольные размеры константы (например, d = & = = 1/// ). Будем решать две задачи Коши с начальными условиями (4) и (5), интегрируя уравнение (2) из 2.3 от поверхности и от дна навстречу друг другу при Z = Л , где 2 т - точка, в которой скорость 0(2) минимальна, то есть С(гт )= min (0(г))
Интегрирование уравнения таким способом, а не напрямую, то есть сразу от поверхности до дна или от дна до поверхности, связано с тем, что полиномы, которыми мы будем аппроксимировать решение, являются функциями, растущими на бесконечности и плохо аппроксимируют экспоненциально убывающие решения. При таком же способе мы всегда аппроксимируем растущие решения, что, как показывает практика вычислений, существенно повышает устойчивость метода. Пусть Пт,н)=і а у (6) ҐІ-0 где Ео может совпадать с 2 , I = 1} . . п , либо лежать в середине отрезка [zc І+І] . Первые два коэффициента Q-n. получаются из начальных условий (4) то есть 0.о = О, Q =6 , а остальные из рекурентных-соотношений а =-; /МгК-е йЛ_з7 (7) полученных подстановкой ряда (6) в уравнение (3). На этом этапе основной вопрос состоит в выборе шага интегрирования, то есть максимального значения [ 2 - о . Для оценки шага интегрирования в уравнении (3) сделаем замену независимого переменного, то есть введем и в результате уравнение (3) примет вид 4 и +(ЩП(АО)+pt) f=0 9) Выпишем связь между производными решения уравнения (9) при Г = О, А07 0:
Произведя вычисление (f+(Yiz) и Ч г ( 2; в точке 2 мы продолжаем интегрирование задачи Коши, снова применяя процедуру разложения в ряд в качестве точки Zo беря 2 и меняя при этом А0 и, если необходимо, то и А1 (при переходе через 2.1 ), так, чтобы вид уравнения (3) не изменился, Начальные данные полу ПТ,г) чаются из вычисленных ПТ: Ч УТ,г) Р путем нормировки р где ИП ЧтЗЧЧ ДЇ V (21) f Yf, г) и % (lj ) вычисля Нормировка начальных данных является существенной особенностью метода, позволяя, с одной стороны, как и в фазовых методах, всегда иметь решения порядка единицы, а с другой стороны, избежать численного интегрирования уравнения Риккати, связанного с большими затратами машинного времени.
Параллельно с вычислением ется функция (22) 0(f) которая является решением неоднородного уравнения (23) - 50 -Для Q+( z) на поверхности Z-0 ставится условие Г(1,о)-о , (1,0) = 0 и на каждом новом шаге интегрирования функция (Ґ(Т Е\ является суммой решений двух начальных задач (24) причем решение первой задачи отыскивается как решение однородного уравнения (3), где вместо If стоит Qf , а решение второй задачи - это решение неоднородного уравнения (23). Аналогично отыскивается У [? ) , /г f? ) а также 0"(;2) » 9iff,z) . причем при Z=H r(\H) fr) $і(ін) о (25) то есть начальные данные для О" есть продифференцированные по f начальные данные для у, В результате при Z - Z т получим условие сшивание решений двух задач Коши (ГЖЫ (?ЛК Ж(?АИ (26) которое является дисперсионным уравнением для рассматриваемой краевой задачи. Для отыскания т мы используем метод Ньютона-Рафсона, то есть L+1 -тая итерация собственного значения получается из I -той с помощью соотношения
Аппроксимация свойств среды (первый способ)...
Для вычисления лучевых траекторий в трехмернонеоднородной среде необходимо задание по измеренным в отдельных точках значениям по крайней мере непрерывной функции трех переменных О (X , U , 2 ) - коэффициенты преломления среды. Если идти по традиционному пути разбиения среды на прямоугольные паралле- лепипеды, внутри которых записывать коэффициент преломления в виде квадратичной формы от трех переменных (X , у , 2 ) [185], то в каждом параллелепипеде необходимо иметь 8 коэффициентов квадратичной формы. Действительно, полная квадратичная форма от трех переменных содержит 10 коэффициентов Р( Ц, )=Рооо+Роо,Х+ Р ноУ+Р оог + Роог + (I) + Pozotf Ргоо?+ Ро»Ху+Р«нХ2 + Р юУг а число условий на р ( л , и , в ) равно 8, то есть числу вершин параллелепипеда, в которых заданы значения коэффициента преломления. Коэффициент преломления, заданный в виде квадратичной формы, обладает, с точки зрения лучевых уравнений, многими удобными свойствами, о которых речь пойдет ниже, однако, прямоугольная форма элементарных блоков, из которых строится среда, вступает в противоречие с произвольной формой дна, а также, если она имеется, с причудливой формой береговой линии океана. Кроме того, для задания значений коэффициента преломления в вершинах прямоугольных параллелепипедов, необходима постановка гидрологических станций в узлах прямоугольной сетки, что представляет определенные трудности в реальных океанических условиях, в то время как почти всегда имеются материалы гидрологических съемок, интересующих нас акваторий, либо хаотически расположенные по акватории, либо, в лучшем случае, вдоль каких либо определенных направлений. В результате учета этих обстоятельств нами было предложено два способа разбиения акватории на элементарные блоки, отличные от параллелепипедов, а именно на вертикальные призмы (настоящий раздел) и симплексы (раздел 3.4), в которых был расписан аппарат лучевой теории.
Коэффициент преломления. Рассмотрим трехмерную среду (Рисунок I). Пусть в точках ІН Ці)) І-1, "М заданы гидрологические станции,то есть массивы значений скоростей звука Си и глубин, на которых они измерены
С помощью описанной выше процедуры триангуляции разобьем множество точек (хс uL) на треугольники и положим эти треугольники в основание вертикальных призм, вдоль ребер которых заданы значения Сі/ . ПОЛОЖИМ Р;;-Пг( с U; 2С/)=тъ? ; А и рассмотрим, например призму с основанием (XiiLJi)l(xlut)Jx5u1 ) Проведем через точки (Xiiyiliij)L-il2lb/j ii.„N1 )il,. ,)l„. 3 плоскости которые разобьют нашу вертикальную призму на меньшие призмы, в каждой из которых зададим Р(х,у,г) в виде то есть в виде квадратичной формы, коэффициенты которой вычисляются при подстановке в (2) значений р(х,Ц,ъ)ъ вершинах призмы: о- р ) vzo р ) ъо- р W - р ) иг- р ) сз р і л Хъ-Хг . X,- Хъ Л-XizlL dl- p -- - ) Ь- р (3) ft = a2(x,, /,,zyj; Pz=nz(Xi,yZlfy); Р Ку ) Непосредственной подстановкой можно убедиться, что форма(2) с коэффициентами (3) принимает заданные значения в вершинах призмы, а на ребрах и гранях призмы определяется лишь значениями в вершинах содержащихся в этих ребрах и гранях и таким образом непрерывно аппроксимирует (х, ,?) при переходе к соседним призмам как по вертикали, так и по горизонтали.
Описанный выше способ аппроксимации коэффициента преломления среды в трехмернонеоднородной среде позволяет получить аналитическую форму лучевых траекторий в пределах одной аппроксимацион-ной призмы и тем самым значительно повысить быстродействие построенных на его основе алгоритмов по сравнению со стандартными методами интегрирования лучевых уравнений, например, как метод Рунге-Кутта.
Поскольку в экспоненты в (21) входит а из (II) имеющее порядок горизонтальных градиентов скорости звука в океане, которые невелики, при программной реализации мошю ограничиться несколькими членами разложения экспонент по степеням аргумента. Если желаемая абсолютная погрешность вычислений есть Є , то, исходя из вида общего члена ряда
Рассмотренный в предыдущем разделе способ аппроксимации свойств среды обладает тем достоинством, что позволяет получить аналитическое решение лучевых уравнений, что сильно ускоряет процедуру вычисления лучевых траекторий. Однако существенным недостатком этого способа является необходимость отыскания корней трансцендентного уравнения при поиске точек пересечения лучевой траектории (формула (21),, раздел 3.3.2.) с границами аппроксимационной призмы. Кроме того поиск корней необходимо производить для всех пяти граней призмы., так как заранее неизвестно через какую из них луч покинет призму. В итоге возникает проблема отыскания такого способа аппроксимации свойств среды, чтобы для нее выполнялись следующие требования: а) непрерывность аппроксимации б) алгебраическая форма лучевых траекторий в) простой анализ условий входа и выхода луча из элементар - 144 -ного блока.
Такую аппроксимацию удалось построить используя некоторые конструкции метода конечных элементов [186] . Она удовлетворяет всем вышеперечисленным требованиям и, кроме того, разбивает среду на кусочно слоистые области, что дает возможность использовать все аналитические соотношения для лучевого метода и метода суммирования гауссовых пучков полученные для слоистой среды в двух- и трехмернонеоднородных средах.
Определение акустической неоднородности среды с помощью звуковых сигналов
В этом разделе описаны особенности, возникающие при численной реализации МСГП. Как было указано выше, на первом этапе из источника выпускается веер опорных лучей, причем на каждом из них производится просчет фундаментальной матрицы ((7), раздел 3.5.1). При переходе луча через разрывы производной С ( 2- ) используется матрица преобразования ((14), раздел 3.5.1) и при отражении от мягкой границы - матрица ((13), раздел 3.5.1). В пределах каждой области непрерывности градиента С ( 2 ) производится анализ наличия каустики по знаку величины Щ г , которая на каустике меняет знак. В итоге, гауссов пучек представляется в виде (слоистая среда): где Qo-Qi+iXQi - 166 К4 - число отражений опорного луча от поверхности, М - число каустик, которых коснулся опорный луч
Подобное представление гауссова пучка позволяет одинаково суммировать вклады пучков, независимо от количества отражений от поверхности соответствующих им опорных лучей и числа каустик, которых эти лучи касаются. Пусть совокупность точек, в которых вычисляется поле равноудалена от источника по горизонтали. Так как каждый гауссов пучек обладает эффективной "шириной , которая в области, где поведение лучей регулярно, равна то локальные координаты, а также величины Q и Г рассчитываются только для точек, которые попадают в "эффективную ширину" гауссова пучка. Таким образом а Д Q - приращение угла скольжения в источнике между соседними лучами, X - расстояние по горизонтали от источника до совокупности точек где вычисляется поле.
Наряду с геометрической характеристикой экспоненциального спада Li , гауссов пучок характеризуется параметром пульсации пучка
Используя эти два параметра можно установить оптимальный режим использования МСГП. Исходя из того, что в регулярной зоне звуковое поле, рассчитанное на основе МСГП, должно сводиться к полю геометрических лучей, мы приходим к тому, что при оптимальном значении параметра JT параметр Li должен достичь своего минимума, а с увеличением частоты МСШ дает верные результаты для более широкого диапазона изменения параметра }f . Минимальное значение "ширины" пучка достигается при Поскольку при расчете поля интеграл заменяется суммой (5), то для снижения погрешности вызванной этой заменой, необходимо достаточное количество опорных лучей, которое должно возрастать с уменьшением значения параметра
Кроме того, число лучей должно быть таким, чтобы сумма (5) отслеживала пульсации величины IX т , характерный размер которых есть L z , поэтому при выборе числа лучей следует руководствоваться минимальным из этих параметров. Отметим, что при увеличении частоты величины параметров падают, что приводит к увеличению надежности МСШ, однако необходимое число опорных лучей возрастает.
При вычислении звукового поля при помощи МСГП следует учитывать тот факт, что на данной частоте излучения могут существовать области, в которых значение Li весьма велико, что есть следствие использованного в разделе 3.5.1 метода стационарной фазы. Это приводит к тому, что в этих областях, несмотря на регулярность поведения лучей интеграл по гауссовым пучкам не может быть сведен к лучевому решению. Аналогичная ситуация может возникнуть когда в одну точку приходит несколько лучей. Хотя в отдельности для каждого из них можно указать приемлемый диапазон изменения параметра гауссова пучка )f , для недостаточно высокой частоты оптимального значения для всех лучей, приходящих в точку не существует. Подобные ситуации устраняются с увеличением частогы излучения. На практике существуют такие расстояния от источника до исследуемой точки, на которых МСГП, ввиду указанных выше причин, неприемлем, однако, если в изменении с расстоянием параметра Li прослеживается некоторая периодичность можно использовать и на значительных расстояниях от источника.
В настоящем разделе мы приведем два примера, иллюстрирующие работу алгоритмов расчета звукового поля лучевым методом в трех-мернонеоднородной среде и МСГП в слоистой среде.
В задаче с применением лучевого метода на первом этапе была построена аппроксимация формы морского дна. Дно задавалось по точкам {tiMi.ii), которые считывались с карты акватории и заносились в память ЭВМ. Массив точек триангулировался с помощью вышеописанной процедуры триангуляции и в каждом треугольнике строилось аппроксимирующее выражение, описанное в разделе 3.2. В точке была построена гладкая поверхность, изолинии которой изображены на рисунке I. На данном рисунке использовался начальный массив из 100 точек.
