Boundary value problem of calculating ray characteristics of ocean waves reflected from coastline

Full Text

Краевая задача о расчете лучевых характеристик океанических волн, отраженных от береговой линии ^[1]

Введение

Ранее в работах [1$–$4] был предложен аналитико-численный подход, осуществляющий быстрый расчет основных параметров длинных океанических волн в рамках теории мелкой воды, состоящий их двух основных этапов: нахождения фронта волны на основе лучей, рассчитанных для заданных граничных условий, и определение высоты профиля с использованием асимптотических формул оператора Маслова. В представленном аналитико-численным подходе, с вычислительной точки зрения, наиболее трудоемким является первый этап $–$ решение краевой задачи. Задача значительно усложняется в случае сложной формы дна океана, а также рассмотрения лучей, отраженных от береговой линии. Наиболее распространенный подход к решению краевой задачи основан на сведении исходной граничной задачи к задачи Коши для гамильтоновой системы лучевых уравнений.

Запишем систему Гамильтона [5, 6] $r\mathrm{<mo>=</mo>}R\mathrm{<mo\; stretchy="false">(</mo>}\varphi \mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo>,}p\mathrm{<mo>=</mo>}P\mathrm{<mo\; stretchy="false">(</mo>}\varphi \mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo>}$ в следующем виде:

$\dot{r}\mathrm{<mo>=</mo>}{H}_p\left(r\mathrm{,}p\right)\mathrm{<mo>=</mo>}c\left(r\right)\frac{p}{\left|p\right|}\mathrm{,}\dot{p}\mathrm{<mo>=</mo>}-H_r\left(r\mathrm{,}p\right)\mathrm{<mo>=</mo>}-\nabla c\left(r\right)\left|p\right|\mathrm{,}$ (1)

где $H\mathrm{<mo\; stretchy="false">(</mo>}r\mathrm{,}p\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo><mo>|</mo>}p\mathrm{<mo>|</mo>}c\mathrm{<mo\; stretchy="false">(</mo>}r\mathrm{<mo\; stretchy="false">)</mo>}$ $–$ гамильтониан, $c\mathrm{<mo\; stretchy="false">(</mo>}r\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}\sqrt{gD\mathrm{<mo\; stretchy="false">(</mo>}r\mathrm{<mo\; stretchy="false">)</mo>}}$ $–$ скорость волны в среде, $g$ $–$ ускорение свободного падения, $D\mathrm{<mo\; stretchy="false">(</mo>}r\mathrm{<mo\; stretchy="false">)</mo>}$ $–$ глубина океана в точке с радиус$–$вектором $r$. Лучи испускаются из начальной точки $r{\mathrm{<mo>|</mo>}}_{t\mathrm{<mo>=</mo>0}}\mathrm{<mo>=</mo>}{r}_0\mathrm{<mo>=</mo><mo\; stretchy="false">(</mo>}{x}_0\mathrm{,}{y}_0\mathrm{<mo\; stretchy="false">)</mo>}$ (точки, где локализовано начальное возмущение ${\eta }^0\mathrm{<mo\; stretchy="false">(</mo><mo\; stretchy="false">(</mo>}r-r_0\mathrm{<mo\; stretchy="false">)</mo><mo>/</mo>}l\mathrm{<mo\; stretchy="false">)</mo>}$ ) с начальными импульсами, лежащими на единичной окружности $p{\mathrm{<mo>|</mo>}}_{t\mathrm{<mo>=</mo>0}}\mathrm{<mo>=</mo>}n\mathrm{<mo\; stretchy="false">(</mo>}\varphi \mathrm{<mo\; stretchy="false">)</mo><mo>=</mo><mo\; stretchy="false">(</mo>}\cos \varphi \mathrm{,}\sin \varphi \mathrm{<mo\; stretchy="false">)</mo>}$.

Запишем краевые условия нахождения траекторий, выпущенных с неизвестным углом из заданной точки $r_0$ и приходящих за некоторое (также неизвестное) время в заданную точку наблюдения $r_1$:

$R\mathrm{<mo\; stretchy="false">(</mo>}\varphi \mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}{r}_1\mathrm{.}$

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

Альтернативным подходом к решению краевой задачи является вариационный расчет лучей на основе принципа Ферма [7]. Ранее в работе [8] был предложен вариант вариационного метода для расчета лучевых траекторий, соответствующих экстремалям функционала времени прихода волны:

$\delta T\left[\gamma \right]\mathrm{<mo>=</mo>}\delta {\displaystyle \underset{A}{\overset{B}{\int }}}\frac{1}{c\left(r\right)}dl\mathrm{<mo>=</mo>0,}$ (2)

где $A$ и $B$ $–$ положения начальной и конечной точек лучевой траектории $\gamma$, $c\mathrm{<mo\; stretchy="false">(</mo>}r\mathrm{<mo\; stretchy="false">)</mo>}$ $–$ скорость волны в точке $r$ на лучевой траектории $\gamma$, $dl$ $–$ элемент длины. Основная идея вариационного метода заключается в прямой численной оптимизации функционала (2), нахождению траекторий, на которых функционал принимает стационарное значение. Главным преимуществом вариационного подхода является автоматическое выполнение граничных условий, так как положения конечных точек зафиксированы в процессе расчета луча. Также на примере модельных сред было показано [8], что вариационный метод позволяет быстро рассчитывать все лучи с заданными граничными условиями и имеет широкие перспективы для дальнейшего применения в задаче распространения океанических волн.

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

Целью данной работы является представление нового варианта вариационного метода для расчета лучей, отраженных от берега, и его применение в модельной среде.

Статья организована следующим образом: в разд. 1 приведено описание нового варианта вариационного метода. В разд. 2 представлены результаты расчетов лучей с различным числом отражений в аналитических моделях океанического дна с линейной и нелинейной геометрией берега. Приведено сравнение лучей, полученных вариационным методом, с решениями системы Гамильтона. Заключение и выводы представлены в разд. 3.

1. ВАРИАЦИОННЫЙ МЕТОД

Рассмотрим модификацию вариационного метода, представленного в работе [8, 11], с учетом отражения лучевой траектории от берега. Функционал времени распространения океанической волны вдоль некоторой кривой $\gamma$ имеет следующий вид:

$T\left[\gamma \right]\mathrm{<mo>=</mo>}{\displaystyle \underset{A}{\overset{B}{\int }}}\frac{1}{c\left(r\right)}dl\mathrm{,}$ (3)

где $A$ и $B$ $–$ положения начальной и конечной точек кривой $\gamma$, $c\mathrm{<mo\; stretchy="false">(</mo>}r\mathrm{<mo\; stretchy="false">)</mo>}$ $–$ скорость волны в точке $r$, $dl$ – элемент длины. Будем считать, что кривая $\gamma$ является начальным приближением луча, отраженного $M$ раз от линии берега на всем протяжении от начальной точки $A$ до конечной точки $B$.

Используя метод прямоугольников, запишем приближенное выражение для функционала:

$T\left[\gamma \right]\approx T\left(r\right)\mathrm{<mo>=</mo>}\frac{1}{2}{\displaystyle \sum _{i\mathrm{<mo>=</mo>1}}^{N-1}}\frac{1}{c\left(r_{i+\mathrm{1<mo>/</mo>2}}\right)}\left|r_{i+1}-r_i\right|\mathrm{.}$ (4)

Кривая $\gamma$ представлена в виде кусочно-линейной кривой со следующими параметрами: $N$ $–$ число узлов (точек) с радиус-векторами $r_i$, $i\mathrm{<mo>=</mo>}\left[\mathrm{1,...,}N\right]$ – номер узла, $c\mathrm{<mo\; stretchy="false">(</mo>}{r}_{i+\mathrm{1<mo>/</mo>2}}\mathrm{<mo\; stretchy="false">)</mo>}$ $–$ значение скорости на интервале $r_{i+\mathrm{1<mo>/</mo>2}}\mathrm{<mo>=</mo>}\frac{1}{2}\left(r_{i+1}+r_i\right)$ между узлами $i$ и $i+1$, $r\mathrm{<mo>=</mo>}\left(r_2\mathrm{,}{r}_3\mathrm{,...,}{r}_{N-1}\right)$ $–$ многомерный вектор, компонентами которого являются вектор $r_i$. Здесь и далее под символами, выделенными полужирным шрифтом, понимаются векторы в многомерном пространстве целевой функции (4). Конечные узлы (точки $A$ и $B$ ) зафиксированы по умолчанию, а положения $N-2$ промежуточных узлов не закреплены на протяжении всего процесса поиска оптимума.

Отметим, что приближенный расчет интеграла (3) может быть осуществлен методом трапеций, Симпсона и другими методами. В данном случае использование метода прямоугольников (4) позволяет избежать проблемы сингулярности подынтегральной функции на берегу ( $c\mathrm{<mo>=</mo>0}$ ), поскольку величина $c\mathrm{<mo\; stretchy="false">(</mo>}{r}_{i+\mathrm{1<mo>/</mo>2}}\mathrm{<mo\; stretchy="false">)</mo>}$ рассчитывается на интервалах $r_{i+\mathrm{1<mo>/</mo>2}}$ кусочно-линейной кривой $\gamma$, на каждом из которых выполняется условие $c\mathrm{<mo>></mo>0}$.

Выделим из числа промежуточных узлов набор точек числом $M$, соответствующих предполагаемым точкам отражения луча от берега. Примем, что начальные положения таких точек должны принадлежать кривой $\epsilon$, описывающей береговую линию. Далее будем учитывать, что смещения всех точек отражения осуществляется только вдоль касательных ${\tau }_j$ к кривой $\epsilon$, где $j\mathrm{<mo>=</mo>}\left[\mathrm{1,...,}M\right]$ $–$ номер точки отражения.

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

Далее последовательно рассмотрим способ оптимизации по переменным $r_i$ и $r_j$. В качестве независимых переменных будем считать положение точек отражения $r_j$.

1.1. Оптимизация промежуточных узлов

Для начала рассмотрим процедуру оптимизации целевой функции (4) по переменным $r_i$ в соответствии с методом обобщенной силы [11].

В случае поиска локального минимума движение каждого $i$ -го узла осуществляется под действием фиктивной силы, основанной на модифицированном градиенте целевой функции $F_{\text{min}}$:

$F_{\text{min}}\mathrm{<mo>=</mo>}-{\left.\nabla T\left(r\right)\right|}_{\perp }+F_{\text{sp}}\mathrm{,}$ (5)

где

${\left.\nabla T\left(r\right)\right|}_{\perp }\mathrm{<mo>=</mo>}\left({\left.{\nabla }_{2}T\left(r\right)\right|}_{\perp }\mathrm{,}{\left.{\nabla }_{3}T\left(r\right)\right|}_{\perp }\mathrm{,}\dots \mathrm{,}{\left.{\nabla }_{N-1}T\left(r\right)\right|}_{\perp }\right)\mathrm{,}$ (6)

$F_{\text{sp}}\mathrm{<mo>=</mo>}\left(F_2^{\text{sp}}\mathrm{,}{F}_3^{\text{sp}}\mathrm{,}\dots \mathrm{,}{F}_{N-1}^{\text{sp}}\right)\mathrm{,}$ (7)

здесь $\left({\nabla }_{i}T\mathrm{<mo\; stretchy="false">(</mo>}r\mathrm{<mo\; stretchy="false">)</mo>}\right){\mathrm{<mo>|</mo>}}_{\perp }\mathrm{<mo>=</mo>}\left(\partial T\mathrm{<mo>/</mo>}\partial r_i\right){\mathrm{<mo>|</mo>}}_{\perp }$ $–$ поперечная к траектории компонента градиента целевой функции вдоль $r_i$, $F_i^{\text{sp}}\mathrm{<mo>=</mo>}\kappa \mathrm{<mo\; stretchy="false">(</mo><mo>|</mo>}{r}_{i+1}-r_i\mathrm{<mo>|</mo>}-\mathrm{<mo>|</mo>}{r}_i-r_{i-1}\mathrm{<mo>|</mo><mo\; stretchy="false">)</mo>}{\tau }_i$ $–$ “фиктивная” сила взаимодействия точек траектории или сила упругости [12], действующая на каждый узел $r_i$, $\kappa$ $–$ константа упругости (в наших расчетах $\kappa \mathrm{<mo>=</mo>0.1}$ ), ${\tau }_i$ $–$ касательная к траектории в точке $r_i$. Введение силы упругого взаимодействия между точками кусочно-линейной кривой позволяет равномерно распределять точки вдоль траектории, избегая областей их компактной локализации, и, также, значительно ускорить сходимость луча [13].

В случае поиска седловой точки, движение каждого $i$ -го узла осуществляется под действием фиктивной силы $F_{\text{saddle}}$:

$\begin{array}{l}{\mathbf{F}}_{\text{saddle}}=\left\{\begin{array}{l}\sum _{\mathrm{\eta }=1}^R\left({\left.\nabla T\left(\mathbf{r}\right)\right|}_{\perp }\cdot {\mathbf{Q}}_{\mathrm{\eta }}\right){\mathbf{Q}}_{\mathrm{\eta }}+{\mathbf{F}}_{\text{sp}}\mathrm{,}\text{\hspace{0.33em}если\hspace{0.33em}}{\lambda }_{\mathrm{\eta }}⩾\text{0}\mathrm{,}\\ -{\left.\nabla T\left(r\right)\right|}_{\perp }+2\sum _{\mathrm{\eta }=1}^R\left({\left.\nabla T\left(r\right)\right|}_{\perp }\cdot {\mathbf{Q}}_{\mathrm{\eta }}\right){\mathbf{Q}}_{\mathrm{\eta }}+{\mathbf{F}}_{\text{sp}}\mathrm{,}\text{\hspace{0.33em}если\hspace{0.33em}}{\lambda }_{\mathrm{\eta }}\text{<0}\mathrm{.}\end{array}\right.\end{array}$ (8)

Здесь $\eta$ $–$ индекс отрицательного собственного значения ${\lambda }_{\eta }$ гессиана $H_{i\mathrm{,}k}\mathrm{<mo>=</mo>}{\partial }^{2}T\mathrm{<mo>/</mo>}\partial r_i\partial r_k$ ( $i$ и $k$ $–$ индекс промежуточного узла) и $Q_{\eta }$ $–$ соответствующий отрицательному собственному значению нормированный собственный вектор, $R$ $–$ число отрицательных собственных значений гессиана.

В результате выражения (5) и (8) составляют основу итерационной процедуры сходимости, в ходе которой положения узлов $r_i$ последовательно смещаются к оптимальной форме. Критерием останова оптимизации является падение модуля максимальной силы $\text{max}\mathrm{<mo>\{</mo>}{F}_i\mathrm{<mo>\}</mo>}$ ниже некоторого малого параметра $\sigma$: $\text{max}\mathrm{<mo>\{</mo>}{F}_i\mathrm{<mo>\}</mo><mo><</mo>}\sigma$. В данной работе используется значение параметра $\sigma {\mathrm{<mo>=</mo>10}}^{-9}$.

1.2. Оптимизация точек отражения

Далее рассмотрим процедуру оптимизации целевой функции (4) по независимым переменным $r_j$. Запишем выражение для частной производной функционала в каждой $j$ -й точке отражения, спроецированной на направление касательной к линии берега ${\tau }_j$:

$\frac{\partial T}{\partial {\tau }_j}\mathrm{<mo>=</mo>}\frac{\partial T}{\partial r_j}-\left(\frac{\partial T}{\partial r_j}\cdot n_j\right)n_j\mathrm{,}$ (9)

где $n_j$ $–$ нормаль к линии берега $\epsilon$ в точке $r_j$. Тогда выражение для силы, действующей на каждую $j$ -ю точку отражения, имеет следующий вид:

${\widehat{F}}_j\mathrm{<mo>=</mo>}\frac{dT}{d{\tau }_j}\mathrm{<mo>=</mo>}\frac{\partial T}{\partial {\tau }_j}+{\displaystyle \sum _{i\mathrm{<mo>=</mo>1}}^N}\frac{d{r}_i}{d{\tau }_j}\frac{\partial T}{\partial r_i}\mathrm{,}$ (10)

где ${\tau }_j$ $–$ касательная к линии берега $\epsilon$, $d{r}_i\mathrm{<mo>/</mo>}d{\tau }_j$ $–$ матрица связи переменных $r_i$ и $r_j$.

Поиск стационарных точек по переменным $r_j$ может быть осуществлен аналогично (5) и (8) заменой $\nabla T\mathrm{<mo\; stretchy="false">(</mo>}r\mathrm{<mo\; stretchy="false">)</mo>}$ на $\widehat{F}\mathrm{<mo>=</mo>}\left({\widehat{F}}_1\mathrm{,}{\widehat{F}}_2\dots \mathrm{,}{\widehat{F}}_M\right)$ :

${\stackrel{}{\dot{\mathbf{F}}}}_{\text{min}}-\stackrel{}{\dot{\mathbf{F}}}$, (11)

$\begin{array}{l}{\widehat{F}}_{\text{saddle}}=\left(\begin{array}{c}\sum _{\eta =1}^R\left(\widehat{F}\cdot {\mathbf{Q}}_{\eta }\right){\mathbf{Q}}_{\eta }\mathrm{,}\text{\hspace{0.33em}если\hspace{0.33em}}{\lambda }_{\eta }⩾\text{0}\mathrm{,}\\ -\widehat{F}+2\sum _{\eta =1}^R\left(\widehat{F}\cdot {\mathbf{Q}}_{\eta }\right){\mathbf{Q}}_{\eta }\mathrm{,}\text{\hspace{0.33em}если\hspace{0.33em}}{\lambda }_{\eta }\text{<0}\mathrm{.}\end{array}\right.\end{array}$ (12)

Выражения (11) и (12) составляют основу сходимости точек отражения $r_j$ к оптимальной форме. Аналогично оптимизации по переменным $r_i$, критерием останова является падение модуля максимальной силы $\text{max}\mathrm{<mo>\{</mo>}{F}_j\mathrm{<mo>\}</mo>}$ ниже некоторого малого параметра $\nu$: $\text{max}\mathrm{<mo>\{</mo>}{F}_j\mathrm{<mo>\}</mo><mo><</mo>}\nu$. В данной работе используется значение параметра $\nu {\mathrm{<mo>=</mo>10}}^{-9}$.

1.3. Начальное приближение луча

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

$q_j\mathrm{<mo>=</mo>}p\left(D\left(r_j\right)-D_0\right)g_j$ (13)

до выполнения условия $D\mathrm{<mo\; stretchy="false">(</mo>}{r}_j\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}{D}_0\mathrm{<mo>=</mo>0}$. Здесь $D\mathrm{<mo\; stretchy="false">(</mo>}{r}_j\mathrm{<mo\; stretchy="false">)</mo>}$ $–$ функция глубины дна, $D_0$ $–$ глубина, на которой необходимо закрепить точку $r_j$, $g_j\mathrm{<mo>=</mo>}\nabla D\mathrm{<mo\; stretchy="false">(</mo>}{r}_j\mathrm{<mo\; stretchy="false">)</mo><mo>/</mo><mo>|</mo><mo>|</mo>}\nabla D\mathrm{<mo\; stretchy="false">(</mo>}{r}_j\mathrm{<mo\; stretchy="false">)</mo><mo>|</mo><mo>|</mo>}$ $–$ единичный вектор вдоль градиента $\nabla D\mathrm{<mo\; stretchy="false">(</mo>}{r}_j\mathrm{<mo\; stretchy="false">)</mo>}$ в точке $r_j$, $p$ $–$ константа, регулирующая скорость градиентного спуска (в нашем случае p =1). Важной особенностью сходимости точек отражения $r_j$, является учет зависимости $r_i\mathrm{<mo\; stretchy="false">(</mo>}{r}_j\mathrm{<mo\; stretchy="false">)</mo>}$ в процедуре оптимизации, в результате чего после каждого смещения $r_j$ требуется коррекция положений всех узлов $r_i$.

1.4. Алгоритм оптимизации луча

Ниже приведем алгоритм оптимизации луча с учетом отражения от берега.

1. Инициализация метода:

(a) задание положений точек $A$ и $B$, числа узлов (точек) $N$ и числа точек отражения $M$;

(б) генерация начального приближения траектории (в данном случае прямой отрезок, соединяющий точки $A$ и $B$ );

(в) фиксация всех $M$ точек отражения на береговой линии (выражение (13)).

2. Расчет силы $F_i$ для промежуточных узлов $r_i$ (выражение (5) для минимума/максимума или (8) для седловых точек). Если выполняется условие $\text{max}\mathrm{<mo>\{</mo>}{F}_i\mathrm{<mo>\}</mo><mo><</mo>}\sigma$, то переход к п. 4. Иначе переход к п. 3.

3. Смещение промежуточных узлов $r_i$ вдоль силы $F_i$. Возврат к п. 2.

4. Расчет силы ${\widehat{F}}_j$ для точек отражения $r_j$ (выражение (11) для минимума/максимума или (12) для седловых точек). Если выполняется условие $\text{max}\mathrm{<mo>\{</mo>}{\widehat{F}}_j\mathrm{<mo>\}</mo><mo><</mo>}\nu$, то переход к п. 6. Иначе переход к п. 5.

5. Смещение точек отражения $r_j$ луча вдоль силы ${\stackrel{}{\dot{\mathbf{F}}}}_j$. Возврат к п. 2.

6. Если выполнены оба условия $\text{max}\mathrm{<mo>\{</mo>}{F}_i\mathrm{<mo>\}</mo><mo><</mo>}\sigma$ и $\text{max}\mathrm{<mo>\{</mo>}{\widehat{F}}_j\mathrm{<mo>\}</mo><mo><</mo>}\nu$, то решение найдено. Иначе переход к п. 2.

2. ЧИСЛЕННЫЕ РАСЧЕТЫ ОТРАЖЕНИЯ ЛУЧА ОТ БЕРЕГА

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

2.1. Плоский берег

Применение вариационного подхода для расчета лучей, отраженных от плоского берега, рассмотрим на примере среды, где функция глубины $D\mathrm{<mo\; stretchy="false">(</mo>}y\mathrm{<mo\; stretchy="false">)</mo>}$ задается следующим выражением [14]:

$D\mathrm{<mo\; stretchy="false">(</mo>}y\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}y\mathrm{.}$ (14)

Скорость распространения волны задается формулой Лагранжа $c\mathrm{<mo>=</mo>}\sqrt{gD\mathrm{<mo\; stretchy="false">(</mo>}y\mathrm{<mo\; stretchy="false">)</mo>}}$, где в качестве упрощения примем $g\mathrm{<mo>=</mo>1}$. Для рассматриваемой модельной задачи воспользуемся аналитическими выражениями для лучевых траекторий из работы [14], которые будем считать эталонными:

$X_1\left(\tau \mathrm{,}\psi \right)\mathrm{<mo>=</mo>}{\left(\frac{\tau }{2}\pm 1\right)}^2\mathrm{,}{P}_1\left(\tau \mathrm{,}\psi \right)\mathrm{<mo>=</mo>}{\left(\frac{\tau }{2}\pm 1\right)}^{-1}\mathrm{,}$ (15)

$\begin{array}{c}Q\left(\tau \mathrm{,}\psi \right)\mathrm{<mo>=</mo>}\frac{\tau }{2}\pm \mathrm{1,}\Theta \left(\tau \mathrm{,}\psi \right)\mathrm{<mo>=</mo>1,}\\ X_2\left(\tau \mathrm{,}\psi \right)\mathrm{<mo>=</mo>}Y\left(\tau \mathrm{,}\psi \right)\mathrm{<mo>=</mo>}{P}_2\left(\tau \mathrm{,}\psi \right)\mathrm{<mo>=</mo>}\Xi \left(\tau \mathrm{,}\psi \right)\mathrm{<mo>=</mo>0,}\end{array}$ (16)

при $\psi \in \mathrm{<mo>\{</mo>0,}\pi \mathrm{<mo>\}</mo>}$ (верхний знак отвечает $\psi \mathrm{<mo>=</mo>0}$, а нижний $-\psi \mathrm{<mo>=</mo>}\pi$ ), а при $\psi \notin \mathrm{<mo>\{</mo>0,}\pi \mathrm{<mo>\}</mo>}$ $–$ функциями

$\begin{array}{c}{X}_1\left(\tau \mathrm{,}\psi \right)\mathrm{<mo>=</mo>}\frac{{{\displaystyle \sin }}^2\varphi }{{{\displaystyle \sin }}^2\psi }\mathrm{,}{P}_1\left(\tau \mathrm{,}\psi \right)\mathrm{<mo>=</mo>}\frac{\sin \psi }{\varphi }\mathrm{,}\\ Q\left(\tau \mathrm{,}\psi \right)\mathrm{<mo>=</mo>}\frac{\varphi }{\sin \psi }\mathrm{,}\Theta \left(\tau \mathrm{,}\psi \right)\mathrm{<mo>=</mo>}{{\displaystyle \cos }}^2\varphi \mathrm{,}\end{array}$ (17)

$\begin{array}{c}{X}_2\left(\tau \mathrm{,}\psi \right)\mathrm{<mo>=</mo>}Y\left(\tau \mathrm{,}\psi \right)\mathrm{<mo>=</mo>}\frac{\tau }{2\sin \psi }+\frac{\sin 2\psi -\sin 2\varphi }{2{{\displaystyle \sin }}^2\psi }\mathrm{,}\\ P_2\left(\tau \mathrm{,}\psi \right)\mathrm{<mo>=</mo>}\Xi \left(\tau \mathrm{,}\psi \right)\mathrm{<mo>=</mo>}\sin \psi \mathrm{,}\end{array}$ (18)

где $\varphi \mathrm{<mo>=</mo>}\psi +\frac{\tau \cdot sin\psi }{2}\mathrm{.}$

Здесь функции, описывающие траекторию в стандартных координатах $\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}p\mathrm{<mo\; stretchy="false">)</mo>}$ в фазовом пространстве Ф обозначены через $x\mathrm{<mo>=</mo>}X\mathrm{<mo\; stretchy="false">(</mo>}\tau \mathrm{,}\psi \mathrm{<mo\; stretchy="false">)</mo>,}p\mathrm{<mo>=</mo>}P\mathrm{<mo\; stretchy="false">(</mo>}\tau \mathrm{,}\psi \mathrm{<mo\; stretchy="false">)</mo>}$, а в окрестности берега в расширенном фазовом пространстве ${\Phi }_{\infty }$ через нестандартные $\mathrm{<mo\; stretchy="false">(</mo>}q\mathrm{,}y\mathrm{,}\eta \mathrm{,}\xi \mathrm{<mo\; stretchy="false">)</mo>}$: $q\mathrm{<mo>=</mo>}Q\mathrm{<mo\; stretchy="false">(</mo>}\tau \mathrm{,}\psi \mathrm{<mo\; stretchy="false">)</mo>,}y\mathrm{<mo>=</mo>}Y\mathrm{<mo\; stretchy="false">(</mo>}\tau \mathrm{,}\psi \mathrm{<mo\; stretchy="false">)</mo>,}\theta \mathrm{<mo>=</mo>}\Theta \mathrm{<mo\; stretchy="false">(</mo>}\tau \mathrm{,}\psi \mathrm{<mo\; stretchy="false">)</mo>,}\xi \mathrm{<mo>=</mo>}\Xi \mathrm{<mo\; stretchy="false">(</mo>}\tau \mathrm{,}\psi \mathrm{<mo\; stretchy="false">)</mo>}$. На фиг. 1 представлен пример семейства лучевых траекторий и фронтов волны, полученных на основе выражений (15) – (18) для набора начальных углов излучений в диапазоне $\mathrm{<mo\; stretchy="false">[</mo>0,2}\pi \mathrm{<mo\; stretchy="false">]</mo>}$ с шагом $\frac{\pi }{10}$ .

 

Фиг. 1. a – Семейство лучевых траекторий (зеленые линии) – решений гамильтоновой системы (1), отраженных от береговой линии; б – фронты океанической волны, полученные на основе лучевых траекторий для моментов времени $T_1$ – $T_6$ с шагом $\Delta T\mathrm{<mo>=</mo>1.2}$. Источник расположен в точке с координатами $\mathrm{<mo\; stretchy="false">(</mo>0,1<mo\; stretchy="false">)</mo>}$. Начальные углы излучений заданы в диапазоне $\mathrm{<mo\; stretchy="false">[</mo>0,2}\pi \mathrm{<mo\; stretchy="false">]</mo>}$ с шагом π/10. Скорость распространения волны $c\mathrm{<mo>=</mo>}\sqrt{gD\mathrm{<mo\; stretchy="false">(</mo>}y\mathrm{<mo\; stretchy="false">)</mo>}}$, где $D\mathrm{<mo\; stretchy="false">(</mo>}y\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}y$. Берег представлен в виде горизонтальной линии серого цвета.

 

Рассмотрим лучевые траектории с заданными граничными условиями $A\mathrm{<mo\; stretchy="false">(</mo>0,<mi>\; </mi>1<mo\; stretchy="false">)</mo>}$ и $B\mathrm{<mo\; stretchy="false">(</mo>5,<mi>\; </mi>1<mo\; stretchy="false">)</mo>}$. В данном случае решениями граничной задачи являются три отраженных луча (см. фиг. 2), соответствующие трем стационарным точкам целевой функции (4). Выполним проверку соответствия полученных лучевых траекторий стационарным точкам целевой функции на основе численной процедуры “экспресс анализа” [15]. Численное исследование функционала (4) основано на понятии “виртуальный луч”. Под термином “виртуальный луч” будем понимать траекторию, разделенную точками отражения на интервалы, где каждый из которых минимизируется независимо. Таким образом, после расчета набора виртуальных лучей для различных положений точки отражения $x_1$ в диапазоне $x_1\mathrm{<mo>=</mo><mo\; stretchy="false">[</mo>0<mo>;</mo>5<mo\; stretchy="false">]</mo>}$ (см. фиг. 2а), была определена функциональная зависимость времени распространения вдоль каждого виртуального луча от положения точки $x_1$ (см. фиг. 2б). Полученная зависимость характеризуется наличием трех экстремумов: локальный минимум 1 (фиг. 2б) и соответствующий ему луч 1 (фиг. 2а), максимум 2 $–$ луч 2, локальный минимум 3 $–$ луч 3.

 

Фиг. 2.a$–$Сравнение лучей 1, 2 и 3, полученных из аналитических выражений [14] (сплошные линии) и рассчитанных вариационным методом (кружки). Серыми тонкими линиями представлен набор виртуальных лучей, используемых в процедуре экспресс-анализа. Берег представлен в виде горизонтальной линии серого цвета; б$–$зависимость времени распространения от положения точки отражения виртуальных лучей. Локальные минимумы и максимум обозначены цифрами 1, 2 и 3.

 

Применение вариационного метода для определения луча 1 осуществлялось на основе прямой минимизации целевой функции. Расчет лучей 2 и 3 вариационным методом осуществлялся в результате последовательного перехода от минимума 1 в минимум 3 через максимум 2 с помощью алгоритма глобальной оптимизации [11]. Полученные результаты согласуются с аналитическими решениями (см. фиг. 2).

Рассмотрим пример расчета лучевых траекторий с несколькими отражениями для граничных условий $A\mathrm{<mo\; stretchy="false">(</mo>0,1<mo\; stretchy="false">)</mo>}$ и $B\mathrm{<mo\; stretchy="false">(</mo>10,1<mo\; stretchy="false">)</mo>}$. Лучевые траектории с одним отражением определены аналогично предыдущему примеру (см. фиг. 2). Лучевые траектории с двумя и тремя отражениями получены в результате последовательного перехода от минимума к максимуму через седловые точки. Результаты расчетов согласуются с решениями гамильтоновой системы (см. фиг. 3). Для случая с двумя отражениями, аналогично предыдущему примеру (см. фиг. 2б), была применена процедура “экспресс анализа”. На фиг. 3 представлена зависимость времени распространения от положения двух предполагаемых точек отражения $x_1$ и $x_2$ вдоль виртуальных лучей. В данном случае задание виртуальных лучей осуществлялось в результате последовательных смещений точек $x_1$ и $x_2$ в диапазонах $x_1\mathrm{<mo>=</mo><mo\; stretchy="false">[</mo>0<mo>;</mo>10<mo\; stretchy="false">]</mo>}$ и $x_2\mathrm{<mo>=</mo><mo\; stretchy="false">[</mo>0<mo>;</mo>10<mo\; stretchy="false">]</mo>}$. Отметим, что в случае совпадения $x_1$ и $x_2$ виртуальный луч вырождался в траекторию с одним отражением. Симметричность полученной зависимости на фиг. 4 связана с особенностью задания виртуальных лучей, а именно возможностью выбора направления движения вдоль луча и возможными перестановками точек отражения $x_1$ и $x_2$. На фиг. 4 представлены экстремумы, соответствующие всем лучевым траекториям: два минимума и максимум, расположенные на оси симметрии, соответствуют лучам с одним отражением (см. фиг. 3а), оставшиеся два минимума, четыре седловые точки и два максимума соответствуют лучам с двумя отражениями (см. фиг. 3б).

 

Фиг. 3.а$–$Лучи с одним, б$–$двумя, в$–$тремя отражениями от берега. Решения, полученные вариационным методом, представлены кружками. Аналитические решения представлены сплошными линиями.

 

Фиг. 4. Карта зависимости времени распространения вдоль виртуальных лучей от положений двух точек отражения $x_1$ и $x_2$. Стационарные точки соответствуют лучам 1 – 7 из фиг. 3a, б.

 

2.2. Круглый бассейн с параболическим дном

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

$D\left(r\right)\mathrm{<mo>=</mo>}a\left({‖r‖}^2-R^2\right)\mathrm{,}$ (19)

где $a\mathrm{<mo>=</mo>}-1$ и $R\mathrm{<mo>=</mo>2}$.

Особенностью применения вариационного метода для данного примера является оптимизация точек отражения $r_j$ вдоль единичного вектора касательной ${\tau }_j$ к берегу, рассчитываемого по формуле

${\tau }_j\mathrm{<mo>=</mo>}\frac{q-\left(q\cdot n_j\right)n_j}{‖q-\left(q\cdot n_j\right)n_j‖}\mathrm{,}$ (20)

где $q$ $–$ произвольно заданный вектор, $n_j\mathrm{<mo>=</mo>}\nabla D\mathrm{<mo\; stretchy="false">(</mo>}{r}_j\mathrm{<mo\; stretchy="false">)</mo>}$ $–$ единичный вектор нормали к линии берега в точке $r_j$.

Результаты расчетов лучевых траекторий и фронтов вариационным методом и методом бихарактеристик представлены на фиг. 5. Вариационный алгоритм последовательно определил четыре луча с одним отражением (см. фиг. 5а) и пять лучей с двумя отражениями (см. фиг. 5б), полностью согласующихся с решениями, полученными методом бихарактеристик. На основе полученных результатов осуществлен расчет величины расходимости лучей в окрестности рассматриваемой точки на основе варьирования начальных углов с последующим построением траекторий методом бихарактеристик и численной оценки якобиана расходимости. В табл. 1 представлены результаты расчетов величин времени прихода волны и расходимости в точке $B$. В табл. 1 номер лучей соответствует нумерации на фиг. 5: луч 1 представляет траекторию без отражения от берега, лучи 2, 3, 4 и 5 – траектории с одним отражением от берега, лучи 6, 7, 8, 9 и 10 – траектории с двумя отражениями соответственно. Расчетные величины времени прихода волны вариационным методом и методом бихарактеристик, представленные в табл. 1, показывают высокую степень совпадения. Результаты расчетов расходимости лучей для рассматриваемого примера (см. таблицу 1) показывают, что наиболее значимый вклад в суммарное поле волны в точке наблюдения дают лучи с одним отражением от берега. Величина расходимости лучей с двумя отражениями значительно больше, в результате чего их вклад в суммарное поле волны в точке наблюдения менее значителен. Отметим, что для луча 9, ввиду близкого расположения каустики, значение якобиана $J\approx 0$, в результате чего расходимость в геометрооптическом приближении не может быть корректно определена.

 

Фиг. 5. Результаты расчетов лучевых траекторий между точками $A$ и $B$ с одной (а) и двумя (б) точками отражения в круглом бассейне с параболическим дном. Решения, полученные вариационным методом, представлены кружками. Лучевые траектории, рассчитанные методом бихарактеристик, представлены сплошными линиями. (в) Фронты волны, восстановленные на основе лучевых расчетов методом бихарактеристик, для моментов времени $T_1$ – $T_4$ с шагом $\Delta T=1.2$. Цветовая шкала соответствует распределению функции глубины водоема $D\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}y\mathrm{<mo\; stretchy="false">)</mo>}$. Берег, где выполняется условие $D\mathrm{<mo>=</mo>0}$, представлен сплошной черной кривой.

 

Таблица 1. Результаты расчета времени прихода $T_1$ и $T_2$ методом бихарактеристик и вариационным методом соответственно и расходимости в точке регистрации $B$ для примера с нелинейным берегом (см фиг. 5).

Номер луча	Время. $T_1$	Время. $T_2$	Расходимость. $\mathrm{1<mo>/</mo>}\sqrt{\mathrm{<mo>|</mo>}J\mathrm{<mo>|</mo>}}$
1	1.21	1.211	0.56
2	2.844	2.844	0.9124
3	2.721	2.722	0.721
4	2.721	2.722	0.721
5	4.526	4.526	0.3471
6	6.084	6.084	0.3645
7	7.728	7.729	0.2763
8	6.661	6.662	0.3766
9	4.2	4.2	–
10	6.084	6.084	0.3645

 

2.3. Обсуждение

Результаты расчетов показывают, что успешное применение вариационного метода основано на соответствии решений гамильтоновой системы типам стационарных точек целевой функции. В данной работе было показано, что тип искомого экстремума зависит от выбранного числа отражений луча от берега. Например, в случае одного отражения, целевая функция с одной независимой переменной имеет два экстремума: минимум и максимум. В случае двух отражений число независимых переменных составляет 2, при этом возможны следующие экстремумы: минимум, седловая точка 1-го порядка, максимум. Аналогично для трех отражений возможны: минимум, седловая точка 1-го порядка, седловая точка 2-го порядка и максимум.

Таким образом, отражение от берега приводит к последовательному вырождению стационарного решения: от минимума через седловые точки различного порядка к максимуму. Число отражений луча, являющееся заранее заданным параметром в вариационном методе, предопределяет возможные типы экстремумов. Наблюдаемая последовательность стационарных точек позволяет использовать алгоритм глобальной оптимизации [11] для поиска всех лучевых траекторий краевой задачи с заданным числом отражений. Отметим, представленные расчеты на аналитических моделях океанического дна позволили

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

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

В то же время применение вариантов лучевых методов $–$ квазиклассического приближения, канонического оператора Маслова, также актуально при решении нелинейной задачи в рамках теории мелкой воды. Именно, в недавней работе [4] было показано, что в случае необрушающихся волн решение нелинейной задачи для уравнений мелкой воды в окрестности пологого берега приближенно эффективно выражается через решения линейной системы мелкой воды и, как следствие, через решения волнового уравнения, но с отсутствующими стандартными краевыми условиями. Способ построения асимптотики решения линейной задачи с достаточно плавным переменным дном также недавно развит в работах [16, 17, 14, 18, 19].

Основным результатом работы является демонстрация возможности применения вариационного подхода для расчета лучей, отраженных от берега. Стоит отметить, что предложенный вариационный метод позволяет рассчитывать многократно отраженные от берега лучи (см. фиг. 3 и фиг. 5), причем число отражений является входным параметром программной реализации метода. Как показали численные расчеты якобиана расходимости лучей (см. табл. 1), с увеличением числа отражений геометрическая расходимость увеличивается. В дальнейшем оценка предельного числа отражений при решении двухточечной задачи может быть сделана на основе якобиана расходимости.

ЗАКЛЮЧЕНИЕ

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

1. рассчитать лучевые траектории, многократно отраженные от берега, с заданными граничными условиями. Сравнение с решениями гамильтоновой системы показало полное согласие;
2. представить соответствие лучей экстремумам различного типа: минимумы, седловые точки, максимумы. Количество и типы экстремумов целевой функции зависит от числа отражений луча от берега;
3. определить все множество лучей с использованием алгоритма глобальной оптимизации в результате последовательных переходов между стационарными точками целевой функции.

Представленный метод имеет перспективы применения в моделях с переменной глубиной океанического дна с нелинейной геометрией береговой линии, в том числе на основе осреднения данных батиметрии, что является предметом дальнейших исследований.

Авторы выражают благодарность Доброхотову С. Ю. и Назайкинскому В. Е. за постановку задачи, внимание к работе и полезные советы. Авторы также выражают благодарность Бессарабу П. Ф. за многолетнее плодотворное сотрудничество, а также Бессарабу Ф. С. и Носиковой В. В. за помощь и полезные советы.

 

[1] Работа по реализации вариационного метода выполнена при финансовой поддержке гранта РНФ $№$ 21-71-30011 (Носиков И. А.). Разработка и развитие метода обобщеной силы и численное моделирование выполнено при финансовой поддержке Проекта Минобрнауки РФ $№$ 075-15-2021-583 (Клименко М. В.).

About the authors

I. A. Nosikov

Demidov Yaroslavl State University

Author for correspondence.
Email: ianosikov@wdizmiran.ru
Russian Federation, ul. Sovetskaya, 14, Yaroslavl, 150003

A. A. Tolchennikov

Insitute of Applied Mechanics, Russian Academy of Sciences

Email: tolchennikovaa@gmail.com
Russian Federation, 101 Vernadskogo Ave., Moscow, 117526

M. V. Klimenko

West Department of the Pushkov Institute of Terrestrial Magnetism, Ionosphere and Radio Wave Propagation; St. Petersburg State University

Email: mvklimenko@wdizmiran.ru
Russian Federation, Pionerskaya str., 61, Kaliningrad, 236035; University Embankment, 7-9, St. Petersburg, 1199034

References

Supplementary files